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1 . SUMMARY 


The objectives of this work are (i) to extend the technique of direct 
numerical simulations to turbulent, chemically reacting flows, (ii) to test 
the validity of the method by comparing computational results with laboratory 
data, and (iii) to use the simulations to gain a better understanding of the 
effects of turbulence on chemical reactions. In particular, we address the 
effects of both the large-scale structure and the smaller-scale turbulence on 
the overall reaction rates, examine the relationship between infinite reaction 
rate and finite reaction rate chemistry, and compare some of the results of 
our calculations with existing theories and laboratory data. 

The scope of work to do this involved the following. First, existing 
computer codes were modified to treat the present problem. Next, extensive 
numerical testing of the computer codes was performed. Finally, in order to 
examine the effects of the mixing layer turbulence, both the large-scale 
structure and the smaller-scale turbulence, on the overall reaction rates, a 
sequence of three problems was computed: (i) reactions on a unidirectional 

(one-dimensional) mixing layer, (ii) reactions on a mixing layer experiencing 
large-scale, two-dimensional vortex rollup, and (iii) reactions on a three- 
dimensional turbulent mixing layer. The simulations of the two-dimensional 
mixing layer with vortex rollup are intended to model the large-scale 
structure in the mixing layer, whereas the three-dimensional simulations 
contain both the large-scale structure and the smaller-scale turbulence. 

The numerical testing involved the comparison of computed results with 
exact solutions for a number of different cases. Both rigid body rotation and 
vortex rollup flow fields were used. The work greatly extended the results of 
earlier work on the advection of a passive scalar on a rigidly rotating flow 
field (the color problem) to include also diffusion, chemical reaction, and 
more complex flows. We have found that high accuracy of the spectral methods 
observed in the advection case is also obtained when these further complica- 
tions are present. Our results indicate that spectral numerical methods may 
prove to be useful in the future both for solving the model equations for 
combustor processes as well as for future studies of chemically reacting 
turbulent flows employing direct numerical simulations. 
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The approach of direct numerical simulations allowed extensive examination 
and interpretation of the reaction process. From the one-dimensional simula- 
tions, we found that we could easily compute finite reaction rates near the 
fast reaction limit, that the results were fairly insensitive to the initial 
conditions (for the class of initial conditions computed), and that the 
results were in reasonable agreement with theoretical predictions. 

For the two-dimensional simulations, we found that in all of the cases 
computed, the vorticity field and the product field approximately coincided. 
From the computation of volume averages, we observed the enhancement of the 
overall reaction rate due to the vortex rollup. It also appeared that the 
merging of vortex cores was a more significant mechanism in increasing the 
overall reaction rate than the straining of the reaction interface. Signifi- 
cant species segregation was apparent, so that, for example, the product of 
the average concentrations was approximately equal to and opposite the 
correlation between the fluctuating concentrations. Some of the higher order 
correlations were also examined. 

In the three-dimensional simulations, the contour plots indicated that the 
vortex rollup takes longer to develop, and the vortices and braids are not as 
distinct as in the two-dimensional case. Also, the vortices that develop are 
not strongly correlated laterally. Both the contour plots and the statistical 
results indicated that the spatial segregation was also not as strong as in 
the two-dimensional case, probably due to the weaker vortex rollup as well as 
the effects of smaller-scale , three-dimensional turbulence. 

Comparisons were made between the simulation results and results using 
similarity theory. Approximately linear growth rates of various computed 
length scales, including the mean velocity half-width, the mean vorticity 
thickness, and the mean product thickness, were obtained and were in agreement 
with the theory. Similarity scaling was found to collapse quite well the 
results for the average reactant concentrations, the rms fluctuating reactant 
concentrations, the concentration correlations, the average product concen- 
trations, and the rms fluctuating product concentrations. 

Some limited comparisons were made with laboratory data. Computed 
profiles that were qualitatively similar to corresponding laboratory profiles 
were obtained for the average reactant concentrations, the rms fluctuating 
reactant concentrations, the average product concentrations, the rms 
fluctuating product concentrations, and the concentration correlations. 
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We have made some comparisons with existing theories. One such theory has 
suggested, in addition to using the equations for the average concentrations, 
including the equations for the concentration fluctuations and correlation. 
These equations can be closed by neglecting certain triple moments when 
compared to certain lower order terms. However, our results indicate that the 
triple moment terms are as important as other terms in the equations, so an 
assumption of this type will probably lead to poor predictions. Another 
proposition is to model the mean reaction term, assuming that it will be 
proportional to the average concentration of the lean species divided by a 
turbulent time scale. We have found that such an assumption will only be 
moderately successful if applied to our case. Finally, it has been suggested 
that the concentration correlation of the reacting species can be estimated in 
terms of that for the nonreacting case, which is much easier to model. 

Although this was proposed mainly for statistically homogeneous flows, we find 
that it is a reasonable approximation for our reacting flow simulations. 
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2. INTRODUCTION 


The need for improved combustor efficiency, reduction in pollutant emis- 
sions, the use of alternative fuels, improved component durability, and better 
combustor outflow characteristics demands an improved understanding of and 
predictive capabilities for the various combustor processes. In a jet engine 
combustor, the geometry, flow rates, locations of fuel injectors, and other 
similar features govern the performance of the combustor. The ability to 
model the effects of these features on combustor performance can greatly aid 
in the design of combustors, reducing development time and leading to more 
optimum designs. Modeling can also minimize the extent of testing and aid in 
the definition of a test program and the interpretation of resulting data. 

Modeling the combustor processes is a very complicated task, involving 
numerical models for the combined chemical, thermodynamic, and aerodynamic 
processes. One of the critical features of the aerodynamics, especially with 
regard to its effects on the chemical reactions, is the highly unsteady tur- 
bulent motion in the combustors. When the time scales of the chemical 
reactions are of the same order as, or much less than, the turbulent mixing 
time, which is often the case in combustors, then the reaction rates are 
controlled by the ability of the turbulence to bring the chemical species 
together. Hence, the overall reaction rates depend more on the turbulent 
mixing time and less on the specific reaction rates (Toor, 1962). 

It is clear that proper treatment of the effects of turbulence is usually 
an essential part of a good combustor model. However, at the present time the 
ability to model the turbulence and its effects on combustors is very limited 
(see, e.g., Mellor and Ferguson, 1980). This is especially true for cases 
where the reaction times are not either very fast or very slow compared to the 
turbulent mixing times (see, e.g., O'Brien, 1931; Libby and Williams, 1981). 

In order to treat the complex, unsteady behavior of turbulence, statis- 
tical methods are generally introduced, and one of two main approaches is 
usually followed: the moment equation method or the probability density 

function method. In the first, the moment equation approach, averages (usually 
time averages) of the relevant physical variables are introduced, and the 
governing chemical/aerodynamic equations are averaged. This leads to a 
closure problem, wherein there are more unknown quantities than equations. 

This problem is resolved by introducing ad hoc assumptions relating various 
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unknown quantities and thus closing the system of equations. (See, e.g. , Libby 
and Williams, 1981, or Donaldson and Varma, 1976, for discussions of this 
approach.) The main difficulty with- this approach is in proposing reasonable 
closure models. This is especially true for the reaction rate terms in the 
equations for the chemical species concentrations. For example, the instan- 
taneous reaction rate r for a binary reaction with an Arrhenius rate is 

r = -RpC C_ exp(-T /T ) 

12 a e 

where is the molar concentration of the species i, T g is the temperature 
and T the activation temperature, R is the rate constant, and p the mixture 
density. If the exponential is expanded in terms of (T /T ), then we see that 

3 6 

the average of r depends on an infinite series of averages of the product of 

C, C„ with inverse powers of T . It is very difficult to truncate the series 
12 e 

and propose reasonable models relating these averages to lower order averages 
needed to close the equations (see, e.g., Borghi, 1974). Furthermore, the 
problem becomes even more intractable for more complex reactions. An alter- 
native is to derive equations for the reaction rate term (Donaldson and Varma, 
1976) and to close the equations at a higher order. However, this leads to 
many more unknowns and equations, and more obscure closure assumptions are 
also necessary. 

The second approach, the probability density function (pdf) approach, 
relies on the equation for the (joint) pdf for the relevant physical variables . 
(See O'Brien, 1981, for a review of this approach.) The closure problem for 
the reaction rate term does not arise. However, closure problems for other 
terms, in particular for the mixing term, do appear, and ad hoc closure assump- 
tions are again necessary. A considerable amount of theoretical work has 
addressed these closure problems (e.g., Curl, 1963; Dopazo, 1976; Pope, 1979; 
Janicka and Kollman, 1979), and an appreciable amount of progress has been 
made. The numerical solution of the resulting equations can become very 
difficult. This is because the number of independent variables in this 
approach, which is equal to the sum of the usual independent variables plus 
the dependent variables in the original problem, can be very large. In order 
to avoid this difficulty, Monte-Carlo methods have been introduced (Pratt, 

1976; Pope, 1979). (See Appendix B for a brief discussion of the relationship 
between the present approach and Monte-Carlo methods.) 
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Hybrid methods have also been used. For example, the pdf approach of 
Janicka and Kollman (1979) has employed moment equations for the turbulence 
kinetic energy and dissipation rates in addition to the equation for the pdf. 
In another approach, assumptions have been made for the form of the 
probability density function, which can be used to derive closure assumptions 
for the moment equations (see, e.g., Jones and Whitelaw, 1982). 

The pdf approach can be very useful in the limit of an infinite reaction 
rate, since the conserved scalar approach can often be employed. In this 
limit the statistics of the reacting species concentrations and of the 
thermodynamic variables can for some cases be found in terms of the pdf of a 
conserved scalar (see, e.g., O'Brien, 1971; Bilger, 1980). 

In order to more adequately model the turbulence and its effects on 
combustor performance using either the moment equation approach or the pdf 
approach, there is a need for better information concerning the turbulence/ 
chemistry interaction. This will allow the formulation of better models and 
more adequate testing of the proposed models. There is also a need for 
improved numerical methods in order to more accurately compute solutions to 
the model equations. To address these needs is the main thrust of the present 
work. 

In the past 5 to 10 years an alternative computational methodology has 
become available for studying transitioning and turbulent flows. Termed 
direct numerical simulations, this approach involves the numerical solution of 
the detailed evolution of the complex turbulent velocity field. Using very 
efficient numerical methods (e.g., pseudospectral methods), we can accurately 
solve the fully nonlinear (possibly low-pass filtered) equations of motion, so 
that closure assumptions are only necessary (if at all) for the smaller-scale 
motions, which have been filtered out. Statistical data are obtained by 
performing spatial, temporal, and/or ensemble averages over the computed flow 
field. 

This procedure is analogous to performing experiments in the laboratory. 
However, it has the advantages that (i) much more statistical information of 
interest can be obtained (since the entire flow field is known at every time 
step), (ii) parameters can be varied easily, (iii) experimental conditions are 
more controllable, and (iv) the effects of large-scale structures can be 
directly addressed. The technique also offers the advantage of circumventing 
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the closure problem. The main disadvantage is the limited spatial and temporal 
resolution available, which limits the range of space and time scales (and 
hence the maximum Reynolds number) that can be included in the simulations. 

At the present time, the approach is usable only as a research tool. 

Direct numerical simulations have been successfully applied to the testing 
of turbulence theories, from analytic theories (e.g., Orszag and Patterson, 
1973; Herring et al . , 1973) to second-order closure modeling assumptions 
(e.g., Herring, 1974; Schumann and Patterson, 1978). These simulations have 
been used to address questions of basic physics of the atmospheric boundary 
layer (Deardorff, 1974), of turbulent diffusion (Riley and Patterson, 1974), 
of mixing layer control through forcing (Riley and Metcalfe, 1980a), and of 
sound generation (Metcalfe and Orszag, 1974). Furthermore, successful valida- 
tion studies have been carried out for homogeneous decay (Mansour et al . , 

1979), turbulent wakes (Riley and Metcalfe, 1980b), mixing layers (Riley and 
Metcalfe, 1980a), and turbulent boundary layers (Moin and Kim, 1982; Patera 
and Orszag, 1981). The approach clearly has potential for use in turbulent 
reaction flows. 

The objectives of the work discussed herein are (i) to extend the 
technique of direct numerical simulations to turbulent, chemically reacting 
flows, (ii) to test the validity of the method by comparing computational 
results with laboratory data, and (iii) to use the simulations to gain a 
better understanding of the effects of turbulence on chemical reactions. In 
particular, we address the effects of both the large-scale structures and the 
smaller-scale turbulence on the overall reaction rates, examine the relation- 
ship between infinite reaction rate and finite reaction rate chemistry, and 
compare some of the results of our calculations with existing theories and 
laboratory data. 

In extending this method to chemically reacting flows, it is important to 
take a step-by-step approach, beginning with simpler problems and then 
proceeding on to more complex ones. We choose as a first step to address the 
problem of an idealized turbulent mixing layer undergoing an irreversible 
binary reaction with no heat release. 

In addressing the objectives of this study, this work will have an impact 
on the main thrusts of the NASA-Lewis Combustion Fundamentals Program as 
enunciated by Mularz (1983), i.e., combustion modeling, model validation, 


7 



fundamental experiments, and numerical methods. With regard to combustion 
modeling, direct numerical simulations is an approach which first, although in 
the near term will only be applied to simplified problems, may in the longer 
term be useful as a combustion model. Second, results of the simulations, in 
analogy to laboratory experiments, can be used for model validation. Third, 
fundamental experiments can be performed with the simulations to elucidate 
basic physical processes. And fourth, as part of the methodology development, 
we have introduced spectral numerical methods into reacting flow problems. 
These numerical methods may prove to be very useful in the future in solving 
the appropriate model equations. 

The scope of work to do this involved the following. First, existing 
computer codes were modified to treat the present problem. Next, extensive 
numerical testing of the computer codes was performed. Finally, in order to 
examine the effects of the mixing layer turbulence, both the large-scale 
structures and the smaller-scale turbulence, on the overall reaction rates, a 
sequence of three problems was computed: (i) reactions on a unidirectional 

(one-dimensional) mixing layer, (ii) reactions on a mixing layer experiencing 
large-scale, two-dimensional vortex rollup, and (iii) reactions on a three- 
dimensional turbulent mixing layer. The simulations of the two-dimensional 
mixing layer with vortex rollup are intended to model the large-scale 
structures in the mixing layer, whereas the three-dimensional simulations 
contain both the large-scale structures and the smaller-scale turbulence. 

In the next section we discuss related work that is especially pertinent 
to the present study. In the third section we explain the methodology of 
applying direct numerical simulations to chemically reacting flows. In the 
fourth section we present our results for the tests of the numerical schemes 
and for the sequence of reaction calculations. Finally, in the last section 
we summarize our results, discussing their implications and also possible 
future directions. 
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BACKGROUND 


Several studies of turbulence using numerical simulations have been carried 
out that are closely related to the present study. Riley and Metcalfe have 
addressed - the development of free turbulent shear flows (without chemical 
reactions) using direct numerical simulations, in particular computing axisym- 
metric turbulent, wakes (Riley and Metcalfe, 1980b) and turbulent mixing layers 
(Riley and Metcalfe, 1980a; Metcalfe and Riley, 1981). Since the present 
study is an extension of this work, it will be briefly reviewed. Also Nielsen 
and Hill (1977; see also Hill, 1979), Ghoniem et al. (1982), and Ashurst and 
Barr (1981) have performed numerical simulations of two-dimensional turbulent- 
like flows with chemical reactions. This work will also be briefly discussed. 

Riley and Metcalfe (1980b) applied direct numerical simulations to the 
downstream development of the turbulent wake of an axisymmetric body. The 
purpose of the calculations was to determine whether the simulations could 
accurately portray the physics of this free turbulent shear flow. The 
approach was to use laboratory data to initialize the flow, to compute the 
downstream development of the flow, and to compare the results with data at 
the appropriate downstream distances. The calculations were fully three- 
dimensional and t ime-de pendent . No modeling was used, so there were no 
parameters to adjust. The results showed good agreement between the 
simulations and laboratory data (to within about 5% for the mean velocity and 
10% for the turbulence intensities) and good agreement with similarity 
theory. In addition, several new features, such as intermittency at the wake 
edge., which were not in the initial conditions, developed in the simulations. 

A similar approach was used to study a turbulent mixing layer (Riley and 
Metcalfe, 1980a), again for the purpose of validation of the methodology. The 
mixing layer calculation was initialized using the laboratory data of Wygnanski 
and Fiedler (1970). The temporally growing layer, rather than the spatially 
growing layer (as in the laboratory experiments), was computed. This was done 
because the numerical methods are much simpler for the temporally growing 
layer, while the physical processes are much the same. Figure 1, taken from 
Riley and Metcalfe (1980a), is a typical sequence of plots of constant 
contours of lateral vorticity, which display the rollup of the layer. This 
compares qualitatively with photographs from laboratory data (see, e.g., 
Chandrsuda et al., 1978). 
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Some idea of the quantitative behavior of the mixing layer simulations can 
be seen by examining the statistical properties of the mean velocity profile. 
Laboratory data and similarity theory show that the temporal growth of the 
mixing layer, as measured by its mean velocity profile half-width z^, should 
be linear. Figure 2 (Metcalfe and Riley, 1981) shows the growth of z M for 
calculations on two different computational meshes (using two different 
computer codes). In both cases the linear growth is exhibited, although at 
slightly different rates. Brown and Roshko (1974) have estimated the spatial 
growth rate of z^ from laboratory experiments. If we assume that space and 
time are related by the transformation 

x = j (U x + U 2 )t 

then these data give 

h 37 2 m ’ °- 0275 

which is close to the values obtained from the simulations. Here is the 
free-stream speed of the high-speed layer, is the free-stream speed of 
the low-speed layer, and U = is the velocity difference across the 

mixing layer . 

The spatial development of the mixing layer can be further studied by 

examining the development of the mean velocity, nondimensionalized by U, as a 

function of transverse distance, nondimensionalized by z^. Similarity 

theory and laboratory experiments show that this scaling should collapse the 

data. Figure 3 (Metcalfe and Riley, 1981) shows that both simulations exhibit 

3 

this behavior. The results for the 64 case show that the data collapse is 
maintained over a time period in which the layer increases by a factor of 
over 5 . 

In addition to addressing the question of validation of the methodology, 
Riley and Metcalfe (1980a) also used the simulations to resolve questions of 
anomalous countergradient momentum fluxes observed in the mixing layer 
experiments of Wygnanski et al. (1980) and Ho and Huang (1982). Using the 
three-dimensional simulations discussed above, and also two-dimensional 
calculations of multiple vortex rollup, Riley and Metcalfe concluded that the 
anomalous results were due to the suppression of subharmonic modes by the 
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harmonic forcing used in the experiments. The experiments by Ho and Huang 
(1982) have supported this conclusion. 

Some work has also been carried out addressing chemical reactions occurring 
on two-dimensional, turbulent-like flows. Nielsen and Hill (1977) computed 
the mass conservation equations for chemical species undergoing convection, 
diffusion, and irreversible chemical reaction on two-dimensional, homogeneous 
turbulent velocity fields. Reasonable agreement was obtained between simula- 
tion results and O'Brien's (1969) independence hypothesis for single species 
reactions. Toor (1962) hypothesized that the decay of the covariance of the 
concentrations for a two-species reaction is insensitive to the rate of 
reaction if the reactants have equal diffusivities and are present in 
stoichiometric proportions. This hypothesis was also successfully tested in 
this study. 

Ghoniem et al. (1982) studied the turbulent combustion in a lean 
propane-air mixture in the mixing layer behind a backward-facing step. They 
used Chorin's (1973) random vortex method to compute the time development of 
the two-dimensional, large-scale structures in the flow field. The flame was 
treated as a constant-pressure deflagration acting at the interface between 
two media and propagating locally at a prescribed normal burning velocity. 

Their results were compared to visualizations from a corresponding laboratory 
experiment using high-speed schlieren photography to observe the flame 
fronts. The comparisons indicated good qualitative agreement between the 
experimental and simulated results. 

Some related simulations were carried out by Ashurst and Barr (1981). 

They also used vortex methods to compute the two-dimensional aspects of a 
mixing layer, but used the flux-corrected transport scheme (Zalesak, 1979) and 
the conserved scalar technique (assuming infinite reaction rates) to treat the 
reaction. Their results demonstrated, in particular, the dramatic effect of 
the large-scale structure in the mixing layer on the diffusion flame 
propagation. 

Some mainly theoretical studies carried out by Marble (1982) and his 
coworkers (see also Karagozian, 1982, and Norton, 1983) are also relevant to 
our work. Using asymptotic analyses for very fast reactions and hence very 
narrow flame fronts, they studied the effects of an isolated vortex on a 
laminar flame front. The isolated vortex can be considered to be a simple 
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model for the large-scale structure in a turbulent mixing layer. They found 

that the augmentation of the fuel consumption due to the vortex is propor- 
2/3 1/3 

tional to T D , and the radius of the core of combustion products is 
proportional to Here, T is the circulation of the 

vortex, and D the reactant diffusivity. 


12 



4. METHODOLOGY 


In order to extend the method of direct numerical simulations of 
turbulence to chemically reacting flows, we start by addressing a turbulent 
mixing layer with the following characteristics. 

i. Binary, single-step, irreversible chemical reaction (A + nB -*■ product) 
with negligible heat release. This implies (a) no temperature 
effects on the reaction rates and (b) no effects of the reactions on 
the flow field. 

ii. Very small Mach number flow so that, together with (i.b), the flow is 
incompressible. 

iii. The statistical properties of the flow are invariant with respect to 
translations in the lateral (y) and flow (x) directions. (See 
Figure 4 for the problem geometry in the calculations.) That is, the 
statistical properties of the flow depend only on the transverse 
coordinate (z ) and time. This enables us to compute temporally 
developing mixing layers instead of spatially developing ones usually 
studied in the laboratory. 

iv. The initial conditions for the reacting species are that the 
reactants are not preraixed, with species A in the bottom half and 
species B in the top half of the flow (see Figure 4). 

The first characteristic allows us to study the effects of turbulence on 
the overall chemical reation rates without the complication of the reaction 
influencing the turbulence. Although the latter element is critical in 
treating combustion problems, it is important in beginning the application of 
direct numerical simulations to chemically reacting flows that a step-by-step 
approach be taken and not too many complications be introduced at one time. 
This is also the reason for the second characteristic. Several laboratory 
experiments carried out in recent years have been subject to the same 
restrictions (i.e., i and ii) and have provided data with which to compare the 
results of our simulations (e.g., Konrad, 1977; Breidenthal, 1978; Mungal , 

1983; and Choudbury et al., 1983). 

The temporally growing layer (characteristic iii) rather than the 
spatially growing layer is treated because the numerical methods are much 
simpler for the temporally growing layer, while the physical processes are 
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much the same. For example, in order to treat the spatially growing layer 
numerically, inflow and outflow boundary conditions are needed, which are very 
difficult to correctly implement for this problem. In the temporally growing 
case, periodic boundary conditions can be used, which are very simple to 
implement, especially using spectral numerical methods. For coflowing mixing 
layers as X = (U^ - + U^) becomes small, where and are the 

free-stream velocities on the high- and low-speed sides of the layer, the 
differences between the two approaches should also become small. Finally, by 
addressing the temporally growing case, the problem can be defined as statisti- 
cally homogeneous in two directions (x and y), so that averages can be taken 
in these two directions, increasing the statistical confidence in the results. 

The above assumptions imply that the (molar) concentrations of the 
reactants, say and C^, satisfy the following diffusion-reaction 
equations: 

it C A + H-TC a - -RC a C b . D a 7 2 C a ( 1) 

It S * S- TC B ■ -" RC A C B * V 2 °b (2) 


Here the stoichiometric coefficient n, the reaction rate coefficient R, and 
the dif fusivities D^ and are all constants. The velocity field u satisfies 
the incompressible Navier-Stokes equations: 


_ 3 _ 

at 


u + u*Vu = 


1 2 

Vp + vV u 

P 


(3) 


V*u = 0 


(4) 


where p is the pressure, p is the constant density, and v is the constant 
kinematic viscosity. 

For the reaction defined above, the molar concentration of the product 
field , Cp, satisfies 


C 

at p 


+ u 


•VC P - rc a c b 


D V 2 C 
P P 


(5) 


where Dp is the molecular diffusivity. 
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For lateral boundary conditions, we define the flow field and reactants to 
be periodic in the flow (x) and lateral (y) directions with periods and 
Ly, respectively. This is consistent with statistical homogeneity in these 
directions. It has been found (see, e.g., Riley and Metcalfe, 1980b) that if 
the length of periodicity is large with respect to integral spatial scales, 
then the influence of the periodicity assumption is minimal. In the 
transverse (z) direction, we define free-slip boundary conditions. In 
particular, these imply that 


3_ 

3z 


u(x,y,z) 


z 



( 6 ) 


37 v(x «y> z) 


= 0 


z = 0, L 


w(x,y ,0) = w(x,y,L z ) = 0 


(7) 

( 8 ) 


Here is the length of the computational domain in the z-direction. 

Free-slip conditions have less influence on the flow field than no-slip 
conditions, since the lateral velocities are not required to be zero at the 
boundaries, and are easy to implement numerically. 

In order to solve the governing equations of motion subject to the 
prescribed boundary conditions, we use pseudospectral numerical methods. We 
choose these methods because in related problems they have been shown to be at 
least twice as accurate in each spatial dimension compared to finite 
difference schemes with the same resolution (Peyret and Taylor, 1983; 
Haidvogel, et al., 1980), and they are competitive in computational 
efficiency. Furthermore, they are very accurate in treating phase information 
(Orszag, 1971), which is of particular importance in computing chemical 
reactions. One of the benefits of this study is the first application of 
spectral numerical methods to chemically reacting flows. 

Pseudospectral numerical methods as applied to this problem involve 
expanding the dependent variables in Fourier-s ine/cosine series, i.e., 


{ u ’ v ’ c a’ c b} = 2 
Ik, I <K 


E 


l |k 2 1 <K 2 0<k 3 <K 3 . 


{ G ’^’ £ A’S} eXp 


/k x k y\ 
2 + 


cosln- 


k 3 z' 


(9) 
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w - 


LEE 

lk ][ l<K 1 lk 2 1 <K 2 0<k 3 <K 3 


w exp 


2tt i 



sin it 



( 10 ) 


(See Orszag and Pao, 1974, for a detailed discussion of a similar application 
of pseudospectral methods.) Here k is the wave number vector, K is the 
truncation wave number vector, and g denotes the transform of g. Note that 
these expansions satisfy the boundary conditions exactly. Solving Equations 
(l) through (5) by pseudospectral methods involves evaluating the spatial 
derivatives in Fourier space but computing the nonlinear terms in physical 
space. Our numerical algorithm uses fast Fourier transforms to evaluate the 
transforms and inverse transforms, and either Adams-Bashforth or leap-frog 
time differencing on the nonlinear terms and Crank-Nicolson (implicit) time 
differencing on the diffusion terms. 

In our series of numerical tests, to be described in the next section, 
which was conducted prior to the numerical simulations, we found that small 
negative concentrations sometimes developed, due to numerical errors. This 
was especially true when steep gradients in concentration were present. These 
negative regions occasionally caused spurious effects when computing the 
reaction rate terms and sometimes led to numerical instabilities. (See McRae 
et al., 1982, for a discussion of this problem and various methods to 
alleviate it.) In order to avoid this problem, at regular intervals during 
the calculations we imposed the following restriction on the reactant 
concentrations : 


C (x , t ) 


(c(x,t) 


l 0 


C _> 0 
C < 0 


( 11 ) 


Here C is the value of the concentration prior to applying the restriction. 
Other, more sophisticated filtering schemes were considered (e.g., the 
nonlinear filter suggested by Forester, 1977), but this scheme was found 
adequate to maintain numerical stability and accuracy and was easy to imple- 
ment. (See Appendix A for a discussion of the types of numerical difficulties 
that arise in computing reactions and our approach to dealing with them.) 

Two computer codes were developed, a two-spatial-dimension plus time code 
(the 2D code) and a three-spatial-dimension plus time code (the 3D code), 
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which employed almost identical numerical schemes. The 2D code was used 
extensively in the numerical testing portion of our study, for the one- and 
two-dimensional simulations, and for debugging the 3D code. The 2D code 
employs either 64x64 or 128x128 point computational grids, with execution 
times of approximately 0.35 and 1.4 sec/ time step, respectively. These times 
are not representative of an optimized code since additional scalar fields 
were computed, as were detailed conservation statistics at every time step. 

The 3D code employs either 32x32x32 or 64x64x64 point computational grids, 
with execution times of approximately 0.7 and 7.0 sec/time step, respectively. 
The 3D code, developed in cooperation with Steven A. Orszag, has somewhat 
better execution times for a given number of grid points since it has been 
fully vectorized and optimized to run on the NASA-Lewis Cray IS computer. 

We have treated both finite-rate and infinite-rate reactions. To compute 
the finite-rate case, we directly solve Equations Cl) through (5), subject to 
the boundary conditions discussed. The reaction rate constant R must be kept 
small enough that we can adequately resolve the reactant fields, which can 
develop steep gradients when R is large. To solve the infinite reaction rate 
case, we use the conserved scalar approach (see, e-g., Toor, 1962; O'Brien, 

1971; Ashurst and Barr, 1981). Consider the scalar defined by 

6 - C A - C B (12) 

[0 is sometimes referred to as a conserved scalar of the Schvab-Zeldovitch type 
(Williams, 1965).] If we assume that the diffusivities of C. and C w are equal, 
then subtracting Equation (2) from Equation (1) shows that 8 is a conserved 
scalar quantity satisfying 

6 + u*ve = DV 2 e (13) 


where D = D^ = D R . If the reaction rate is infinite, then the two reactants 
cannot coexist at the same physical point, i.e., the species are segregated. 
In this case, the reactants are directly related to the conserved scalar 0 by 


C 


A 


C 


B 



0 > 0 

0 < 0 

0 _> 0 

0 < 0 


(14) 


(15) 
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Therefore, in the fast reaction limit, by computing the time history of the 
conserved scalar 6, we can obtain the entire time history of the reactants. 

We can also obtain the product field in the infinite rate case by 
considering the conserved scalar defined by 

4> - c p + C A (16) 

Assuming that the diffusivity of the product is the same as for both of the 
reactants, then by adding Equations (l) and (5) we find that <j> also 
satisfies Equation (13), i.e., it is a conserved scalar. (Note, however, that 
the initial conditions for <j) are different than those for 0.) Once C . 
has been obtained using 9, then Cp can be obtained from <J> . 

For initial conditions for the reactant concentrations, we use the 
following functional forms: 

z 

c a ( 2 ’ 0) = 1 f ex P "C 2 /z 2 d£ (17) 

z 0 /* 

z 

C fi (x,0) = — f exp -C 2 /z 2 d£ (18) 

z oA 4, 

Note that there are no initial concentration fluctuations. The initial product 
concentration, Cp(x,0), is taken to be identically zero. We choose these 
functional forms because they give profiles not unlike those measured in the 
laboratory (see, e.g., Konrad, 1977), they are smooth and easily resolvable 
using sine/cosine expansions, and analytic solutions can be obtained for 
special cases using these initial conditions. 

In order to separate the effects of the large-scale, quasi-two-dimensional 
structure and the smaller-scale three-dimensional turbulence in a mixing layer 
on the overall reaction rates, we have considered a series of three simula- 
tions. The simulations all have the same initial conditions for concentra- 
tion, but differ in the initial conditions on the velocity fields. The three 
components in the series are (i) laminar, unidirectional flow, (ii) two- 
dimensional vortex rollup, and (iii) three-dimensional flow with vortex rollup 
and smaller-scale turbulence. The two-dimensional simulations model the 
large-scale structures in the mixing layer, and the three-dimensional 
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simulations contain both the large-scale structures as well as the smaller- 
scale turbulence. The initial conditions for these three flows are defined as 
follows . 

For the laminar, unidirectional flow, the initial velocity is taken to be 

u(x,0) = [u(z),0,0 J (19) 


U(z) 


2 tanh 



( 20 ) 


Here U is the velocity difference across the layer, and z^ is the distance 
from the plane of symmetry to the transverse plane where the mean velocity 
rises to one-half its free-stream value. This profile is a good approximation 
to the mean velocity profiles measured in the laboratory in a mixing layer 
downstream of a splitter plate. The time development of this field is very 
simple, since the nonlinear terms in the momentum equation [Equation (3)] are 
identically zero. The profile broadens with time as the velocity diffuses in 
the transverse direction. With the kinematic viscosity set to zero, the 
velocity field remains stationary in time. The chemistry problem is also 
greatly simplified for this case, since the convection terms in the species 
conservation equations are also identically zero. Therefore, mathematically 
the problem reduces to the one-dimensional problem of the diffusion of one 
species into another, which has been studied extensively and for which exact 
solutions are available (e.g., Burke and Schumann, 1928). 

For the two-dimensional vortex rollup case, we define the initial 
conditions for the velocity field as 


u (x,0) 


U(z) “ |r 'Kx), 0 


( 21 ) 


The mean velocity is the same as in the previous case. However, a perturbation 
has been added to it in terms of a stream function \[i defined by 


ijj(x) = ijj^(x) + i[» SH (x) 


( 22 ) 


Here 4'^ is the stream function for the most unstable mode of the mean velocity 
profile, and is the stream function of the subharmonic of the most unstable 

mode. (See Michalke, 1964, for a discussion of the properties of these modes 
and Riley and Metcalfe, 1980a, for results of calculations using these initial 
conditions.) The presence of only the fundamental mode in the mixing layer 
produces a single vortex rollup. When the subharmonic is added in, but 
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out-of-phase* with the fundamental, then a second rollup occurs (see the 
figures discussed in Section 6.2 of this report and Peyret and Taylor, 1983, 
p. 4, for contour plots of the vorticity field during rollup). The 
calculations reported here have used various initial conditions consisting of 
the fundamental mode alone, the subharmonic mode alone, and the fundamental 
and subharmonic modes added together out-of-phase. 

The methodology for initial conditions for the three-dimensional case is 
similar to that used by Riley and Metcalfe (1980a) in their three-dimensional 
simulations. (See also Riley and Metcalfe, 1978, and Orszag and Pao, 1974, 
for a discussion of this initialization procedure.) The initial velocity 
field is written as 

u(x,0) = ^U(z) + u'(x), v'(x), w'(x)J (23) 

where again the mean velocity is taken to be the same as in the two previously 
discussed cases. The fluctuating velocity field u' is defined by 

u'(x) = P jl(z) A(x) | (24) 


where the form function l(z) is used to define the transverse profile of the 
root-mean-square (rms) velocity, and the operator P is a projection operator 
(in Fourier space) used to insure that the initial velocity field is incompres- 
sible. The vector A is written in terms of its Fourier-sine/cosine transform 
[See Equations (9) and (10)], and its spectral amplitudes are chosen so that 
the velocity field u' has the longitudinal (x) energy spectrum 


2 ~2 

E (k., ) = £ u^ 
11 TT 


+ kt A 


(23) 


where A is the longitudinal integral scale and u the turbulence intensity. The 
laboratory data of Wygnanski and Fiedler (1970) were used to initialize the 


*The term "out-of-phase" means that the cores of the fundamental mode vortices 
do not coincide with the cores of the subharmonic vortex. 
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flow. The initial rms profile, which defines the function l(z), was taken to 
be 


u(x) 

■ — ~ = 0. 18 exp 


^-0.147 



(26) 


where z M is as defined above. (Note that the exponential coefficient is 
given incorrectly in Riley and Metcalfe, 1980a, and in Peyret and Taylor, 
1983, p. 257). The longitudinal integral scale A, which defines the energy 
spectrum [Equation (25)] , was taken to be 


— = 2.0 (27) 

Z M 


In the three-dimensional simulations, in order to obtain a single vortex 
rollup, the size of the computational domain in the streamwise direction, 

L x> was taken to be equal to the wavelength of the fundamental mode. To 
obtain a double rollup, L x was extended to be the wavelength of the 
subharmonic . 

In a particular realization of the flow field, the vector A is generated 
numerically, consistent with the energy spectrum defined by Equation (25), 
using a pseudorandom-number generator. The resulting velocity field is a very 
complex, three-dimensional field which, when averaged, has approximately the 
desired mean velocity, turbulence intensity, and longitudinal energy spectrum. 
As previously mentioned, the temporally growing mixing layer is consistent 
with statistical homogeneity in the x-y plane (or along the x-axis in the 
two-dimensional case). Therefore, to obtain statistical quantities from our 
simulations, averages were taken in x-y planes at a given z, i.e.. 


f(x) 


_1 

N N 
x y 


N , N 
x y 



j=l 


f ( iAx , jAy , z ) 


(28) 


Here N x and are the number of grid points in the x- and y-directions , 
respectively, and (Ax,Ay,Az) define the computational grid spacing. These 
averages correspond to averages computed from data obtained from probes fixed 
in space in a laboratory experiment in which the flow moved past the probes. 

In the results presented, velocities are nondimensionalized by U, the mean 
velocity difference across the layer. Time is nondimensionalized by the 
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inverse of the maximum initial mean vorticity, T = (dU/dz) \ and distances are 
nondimens ionalized by the length, L = UT. Species concentrations are normalized 
by = C^, t * ie ^ ree-streain value of C^. Important nondimensional numbers 
for our problem are the Reynolds number, R^ = UL/v, the Schmidt number. Sc = v/D, 
and the first and second Damkohler numbers, D T = LRC /U and D TT = L^RC / D. The 

I oo II CO 

latter describe the ratio of the convection time scale to the reaction time 
scale and the ratio of the diffusion time scale to the reaction time scale. 

Table 1 contains a list of these parameters for the calculations discussed in 
this report. 

Note that, for the cases reported herein, we have, without loss of genera- 
lity, taken the stoichiometric coefficient to be 1. We have only considered 
the restricted case where the diffusivities of both of the reactants and the 
product are the same and have only treated reacting species whose free-stream 
concentrations are in stoichiometric proportion, i.e., = C^. Realizing 

that the numerical resolution requirements for the velocity field and species 
concentration fields are very similar, we have also selected the Schmidt number 
to be 0(1). Note that this is approximately the case in most gaseous mixtures. 
The entire range of reaction rates is considered, from no reaction, through 
moderate speed reactions, to the infinite reaction rate limit. 

The Reynolds number was kept low enough to adequately resolve the fluid 
motion on the computational grid. The maximum turbulent Reynolds number, 
where the turbulent Reynolds number is defined by 



with X the Taylor microscale, was about 50. This is considerably lower than 
turbulent Reynolds numbers for related laboratory experiments. For example, 
in the experiments of Breidenthal (1978), who worked with chemical reactions 
in water, the turbulent Reynolds number was typically the order of several 
hundred. This difference in Reynolds numbers does not appear to be too 
important for quantities that depend principally on the large-scale structure 
of the turbulent flow, e.g., the mean velocity and turbulence intensities (see 

v. 

Riley and Metcalfe, 1980a, and Metcalfe and Riley, 1981). However, for 
reactions occurring on turbulent flows, which depend critically on the small- 
scale turbulence stretching out the reaction zone to enhance the overall 
reaction rate, this difference in Reynolds numbers could be significant. 
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5. NUMERICAL TESTS 


In order to determine the effectiveness of applying pseudospectral 
numerical methods to reacting flow problems and also to validate the computer 
codes, we defined a series of test problems of increasing complexity. The 
first several were based upon the color problem, i.e., the advection of a 
passive scalar field by a velocity field consisting of rigid body rotation. 
Exact, analytical, t ime -dependent solutions for the concentration fields were 
obtained in specific cases. Other tests were carried out for more complex 
flow fields by examining the conserved scalar, i.e., the difference between 
the two computed concentrations of reacting species [see Equation (12)]. The 
velocity field computations were checked by comparing computed solutions with 
solutions obtained from other validated computer codes. Among the critical 
factors affecting the numerical accuracy in computing reacting flows are the 
tendency to develop steep concentration gradients for large reaction rates and 
the spurious effects caused by the presence of negative concentration values, 
which occur near regions containing steep gradients. As mentioned in the 
previous section, a simple filtering scheme [see Equation ( 1 1 ) ] was used to 
avoid this latter problem. 

The first case computed was the color problem (see, e.g., Orszag, 1971; 

McRae et al., 1982). In this problem, initially conically distributed 
concentration fields are rigidly advected in a circular motion around an axis 
of rotation. There is no diffusion or chemical reaction. The exact solution 
at time t is just the initial field rotated through the angle fit, where fi is the 
angular velocity of advection. This is a nontrivial test of the numerical 
scheme, since the equations are solved in a fixed frame of reference so that 
the concentration field is rotated through the computational grid. Figure 5 
is a plot of the initial concentration field for these tests, showing a 
localized, conical distribution. All of the tests using the rigidly rotating 
flow field were performed on a square 64x64 point computational grid with each 
side of length tt, the rotation rate fi was 2tt , and the time stepping was selected 
to give 2000 time steps per revolution so that time-stepping errors were 
negligible. In these first color problem tests, the filtering scheme was not 
used. Figure 6 shows the concentration field after one revolution. Since the 
exact solution after one revolution should just be equal to the initial field, 
a comparison of Figures 5 and 6 shows the distortion of the solution due to 
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numerical errors. Clearly the sharp peak and other features of the cone are 
accurately maintained. The maximum deviation from the exact solution 
(L^-error) is 2.36%. These results are similar to those of Orszag (1971) 
and should be compared with the results of McRae et al. (1982), which were 
obtained using other numerical methods and are of significantly poorer quality. 

To determine the accuracy of computing the reaction rate terms in 
Equations (l) and (2), we computed the rigid rotation color problem just 
described, but with the reaction rate terms active. We chose the initial 
concentration fields to be identical. The exact solution for the 
concentration of either species can be easily found to be 


C(x' ,t) 


C(x'.,0) 

1 + R C(x' ,0) t 


(30) 


where C(x,0) is the initial concentration field, and x' is the position x 
measured in the rotating coordinate system. The peak value of the initial 
concentration was 1.0. Figure 7 is a plot of the computed concentration field 
after one-half revolution with no reaction, while Figure 8 gives the concen- 
tration of either of the reactants at the same time for a case with the 
reaction rate coefficient R equal to 1.0. Due to the reaction process, the 
shape has changed somewhat, and the peak value has dropped to 0.66. Again 
small oscillations are observed throughout the field, with the L^-error equal 
to 2.76% of the initial peak value. Figure 9 shows the results of a similar 
calculation at one-half revolution with R equal to 4.0. The peak has dropped 
to a value of 0.334, and the L -error is 2.64% of the initial maximum. 

CO 

The accuracy of the calculation of the diffusion terms (without reaction) 
in Equations (1) and (2) can be determined by computing the evolution on the 
rigidly rotating flow field of the scalar field with the initial condition 

/ 2 2 

I X 2 

C(x,z,0) = exp I 2 2 
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In the rotating reference frame (X 1 ), the analytical solution is 
C(x',z',t) = — - — - ■ — — exp 
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Figure 10 is a contour plot and Figure 11 is a perspective plot of the initial 
concentration field, and Figures 12 and 13 show the resulting computed 
concentration field after one revolution. The comparison with the exact 
solution shows that the errors are within 1% of the peak value. 

To test the numerical accuracy with both the diffusion and reaction terms 
present, we chose to address the evolution of two species diffusing and 
reacting across a mixing layer and use solutions which could be obtained for 
the conserved scalar [see Equation (12)]. In these and the following 
simulations with chemical reactions, the filtering scheme [Equation (11)] was 
applied every N time steps. In an accurate calculation, the solution should 
not depend strongly on N, and we found that choosing N between 5 and 15 worked 
well. For the initial concentration profiles given by Equations (17) and 
(18), and for the laminar, unidirectional flow field [Equation (19)], the 
exact solution for the conserved scalar field 0 can be found to be 



Using the 2D code with a 64x64 point computational grid, the one-dimensional 
mixing layer problem was computed for a range of dif fusivities and reaction 
rate coefficients (see Table 1). The conserved scalar 0 was computed from 
the difference between the two computed concentration fields and was then 
compared to the exact solution. Errors were less than 1% of the maximum 
concentrations . 

We next applied the same initial concentration profiles, but used the 
two-dimensional vortex rollup velocity field [see Equations (21) and (22)]. A 
64x64 point computational grid was again used. The diffusion reaction 
equations were solved for and Cg on this flow field, and the conserved 
scalar computed from these solutions. The conserved scalar was also computed 
directly by solving the diffusion equation (without reaction). Because of the 
fast convergence of the spectral series for 0, its numerical solution can be 
obtained very accurately when it is solved for directly. The solution for 0 
is independent of the reaction rate R, and so is equally valid when steep 
gradients in and Cg are present. Thus, a comparison of solutions for 
0 computed by the two different methods serves as a useful check of the 
methodology. For run MR04A, Figure 14 gives a contour plot of the conserved 
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scalar computed indirectly from the difference between the species 
concentrations, and Figure 15 gives a similar plot for the conserved scalar 
computed directly, both at the time of 12, at which vortex rollup has 
occurred. Except in the very smallest details of the +1 and -1 contours, the 
plots are virtually identical. These results are typical of the cases 
computed. This approach of using the conserved scalar as a numerical check 
was used throughout our subsequent calculations, both with the 2D and 3D codes. 

We also checked the accuracy of each calculation by computing the species 
conservation invariant 


I = 


/ (v ; S - 2 s) 


dv 


(34) 


In the simulations listed in Table 1, over the course of a calculation this 
quantity generally deviated from the initial value by less than 0.02%. 

The accuracy of the velocity field computation was checked by comparing a 
calculation of nonlinear vortex pairing in a mixing layer with previous 
calculations performed using a different numerical code that had been fully 
debugged and validated (Riley and Metcalfe, 1980a). The results of the 
calculations using the two codes were essentially identical. This completed 
the testing and debugging of the 2D code. 

The 3D code, which was developed in collaboration with Steven A. Orszag, 
was validated by using it to compute two-dimensional problems and comparing 
the results with those from the 2D code simulations. Test runs were made for 
the two-dimensional reacting mixing layer, both with and without vortex 
rollup. In addition, in order to test the three-dimensional capabilities of 
the code, separate calculations were made with the direction of the mean flow 
parallel to each of the horizontal axes. Finally, the velocity field was 
checked by performing simulations of fully turbulent mixing layers and 
comparing the results with previous simulations made with different computer 
codes (Riley and Metcalfe, 1980a; Metcalfe and Riley, 1981). 
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6 . NUMERICAL RESULTS 


6.1 One-Dimensional Mixing Layer 

The first in our series of simulations was the calculation of the time 
development of the concentration of reacting species on a one-dimensional 
mixing layer. The purposes of these simulations were (i) to examine the 
dependence of the overall reaction rate on the reaction rate coefficient R, on 
the diffusivities [see Equations (1) and (2)], and on the initial concentration 
profiles and (ii) to establish a base against which to compare the effects of 
vortex rollup and three-dimensional turbulence on the overall reaction rates. 

The initial conditions for the concentration of the reacting species are 
given by Equations (17) and_(l8), and those for the velocity field are given 
by Equations (19) and (20). As discussed in Section 4, for this case the 
convective terms in the species concentration equations [Equations (1) and 
(2)] and the Navier-Stokes equations [Equation (3)] are identically zero, so 
that the velocity field does not influence the behavior of the concentration 
fields. The problem reduces to that of the interdiffusion and reaction of the 
two species, which has received much attention in the past. Burke and 
Schumann (1928) obtained the exact solution to this problem for the special 
case of an infinite reaction rate and initial concentrations defined by 

C.(z,0) = C. H(-z) (35) 

A Ac° 


C B (z,0) = C Bc H(z) (36) 

where H(z) is the Heaviside step function, and C. and CL are the free-stream 

Aon Bco 

values of and Cg, respectively. For the case of the reactant concentrations 
in stoichiometric proportion away from the reaction zone (i.e., C Aoo = C Boo ), 
their solution for the total (integrated) product becomes 

OC 1/2 

> / C P ( 5' C) d 5 - 2 C a»( ! t) <37) 

“00 

The notation f is used to denote the result obtained by integrating f over the 
entire domain in z. 
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6.2 Two-Dimensional Mixing Layer 

The second in our series of simulations was the two-dimensional mixing 
layer. In these simulations, a small perturbation is added to the initial 
velocity field [see Equations (21) and (22)], which causes the flow field to 
roll up. This vortex rollup is qualitatively similar to that observed in 
laboratory visualizations of laminar and turbulent mixing layers and has some 
quantitative features in common 0 with them as well (Riley and Metcalfe, 

1980a). This two-dimensional vortex rollup appears to be the essential 
large-scale structure in a three-dimensional, turbulent mixing layer (see, 
e.g., Browand and Ho, 1983). The two-dimensional simulations thus allow us to 
separate out the effects of the large-scale structure from the smaller-scale 
turbulence and, at the same time, to take the first step in simulating the 
turbulent mixing layer. The purposes of this step in our series of 
simulations are (i) to isolate these large-scale structures in the mixing 
layer and (ii) to examine their influence on the chemical reactions. 

For these simulations, the same initial conditions for the concentration 
fields are used as for the one-dimensional case [Equations (17) and (18)]. 

The perturbation added to the initial velocity field was in the form of the 
most unstable (fundamental) mode and/or the subharmonic of this mode. The 
mathematical solutions for these modes are obtained by solving the 
Orr-Sommerfeld equation for the velocity profile, given by Equation (20), and 
are discussed in detail by Michalke (1964). 

Figure 19 gives a sequence of constant contour plots of the vorticity for 
the case where the initial perturbation consisted of the most unstable mode 
alone (Case 1). Because vorticity is conserved following a fluid element 
(except for diffusive effects) for two-dimensional flows, vorticity is useful 
as a tracer to visualize the flow. Vortex rollup is discernible at a time of 
4, the rollup attains its maximum extent at about a time of 12, after which 
the vortex collapses since there is no subharmonic present for pairing to take 
place. When the subharmonic alone is present (Case 2, see Figure 20), vortex 
rollup is underway at a time of 12, the maximum extent of rollup occurs at 
about a time of 24, and subsequently the layer collapses. The behavior of 
these two flows is very similar, the characteristic length and time scales 
being a factor of about 2 greater for the subharmonic case (Acton, 1976). 
Finally, when the fundamental and subharmonic are added together out-of-phase 


30 



(Case 3, see Figure 21), a double rollup occurs, the primary rollup occurring 
at about a time of 8, and the secondary rollup peaking at about a time of about 
24. The growth of the mixing layer, as measured by z^, the distance from the 
axis of symmetry to the point at which the average value has attained one-half 
of its free-stream value, is given for each of these cases in Figure 22. The 
general qualitative conclusions obtained from the contour plots are evident. 
(See Riley and Metcalfe, 1980a, for a more detailed discussion of the mechanics 
of the growth and collapse of the vortex cores.) 

The contour plots of the concentrations of the reactant and product fields 
give considerable insight into the reaction history. Figure 23 contains con- 
tour plots of species A and B for Case 1 at the time of 12, when the rollup has 
been completed and collapse has just begun. We see that large regions of fluid 
containing species A have been convected from the upper into the lower region 
and vice versa. The reaction zone is defined by the interface between the two 
species and has been lengthened significantly by the vortex rollup. The 
gradients of the two species appear very steep at this interface. Constant con 
tours of the product field are given in Figure 24 for the same time. We see 
that the product field is located along the reaction zone, that it has very 
steep gradients along the braids of the vortices, and that it is more diffuse 
in the vortex core. It is interesting to compare the constant contour plots of 
the product field and vorticity field (Figure 19d). Both only exist near the 
reaction zone, gradually diffusing outward, so that their appearances are very 
similar. A contour plot of the product field at a time of 24 (Figure 25) shows 
that, after collapse has occurred, mixing and reaction continue in the central 
layer, causing a somewhat uniform concentration of product there. 

Figure 26 contains a sequence cf constant contour plots of the con- 
centration of species A for Case 2, the subharmonic rollup case. The severe 
stretching of the reaction zone and vortex rollup are clearly visible by a . 
time of 24. Again, we see a large region of fluid containing species A 
protruding into that containing species B. At the time of 36, collapse has 
occurred, and species A has almost been depleted in the interior of the vortex 
core. Therefore, although the reaction interface has been significantly 
lengthened by the vortex rollup, it has also been shortened due to the 
depletion of species A (and species B) in the vortex interior. This latter 
effect has been conjectured in several theories (see, e.g., the "flame 
shortening" mechanism discussed by Marble and Broadwell , 1977). Figure 27 is 
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a constant contour plot of the product at the time of 24. The product again 
appears somewhat diffuse in the vortex core but localized along a thin region 
in#the braids. Comparing Figure 27 with the contour plot for the vorticity at 
the same time (Figure 20b), the similarity in the profiles for these 
quantities is again noted. 

Figure 28 contains a sequence of constant contour plots of the concen- 
tration of species A for Case 3, the double rollup case. The rollup of the 
fundamental mode and the initial stages of the subharmonic are evident at the 
time of 12. At a time of 24, the peak in the second rollup has occurred. 
However, species A has been significantly depleted in the vortex core at this 
time, another example of flame shortening. By a time of 36, collapse is 
occurring, and the continued stretching of the reaction interface and flame 
shortening are clearly evident. The fact that flame shortening has happened 
somewhat sooner for Case 3 than Case 2 (compare Figures 26b and 28b) indicates 
that more reactant is being consumed in Case 3, so that the double rollup 
causes a larger overall reaction rate. Figure 29 shows the product field at a 
time of 24. The vortex core, in which the product is rather diffuse, is 
clearly evident, as well as the thin layer of product along the braids. 
Comparing Figure 29 with the corresponding plot for the vorticity (Figure 21b) 
again shows the similarity between the two quantities. 

As in the one-dimensional case, calculations were performed for different 
initial concentration profiles and different reaction rate coefficients (see 
Table l). For the range of concentration profile length scales [z q defined 
in Equations (17) and (18)] studied, the results were again relatively 
insensitive to this length scale (see Figure 18). Furthermore, for the 
reaction rate coefficients considered, the results were close to the infinite 
reaction rate limit, and hence the reactions were diffusion and convection 
limited. 

In order to compare the two-dimensional results with the one-dimensional 
results, integrated (total) averages were defined by 


C(t) 



C(£ , t) d£ 


(38) 


where g defines the integral over all z of the function g, f indicates a 
spatial average of f [see Equation (28)], and the transverse (z) boundaries in 
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the computational domain are denoted by 0 and L^. Note that of course for the 
one-dimensional case, f = f. Figure 30 gives the time development of the inte- 
grated mean product both for the one-dimensional simulations and for Case 1 
(the fundamental mode alone) of the two-dimensional simulations, and for three 
different diffusivities. Comparing a one- and a two-dimensional case for a 
given diffusivity, we see that for short times the total products are about 
the same. However, as the mixing layer rolls up in the two-dimensional 
simulation (at a time of about 8), the total product is greatly increased. 

The final value of the product for these two-dimensional simulations is also 
plotted in Figure 17. The enhanced amount of product due to the vortex rollup 
is evident. Also, the dependence of the results on the diffusivity is 
significantly different. The velocity field has to some extent exerted a 
controlling influence on the overall reaction. 

Figure 31 gives a similar plot for the integrated mean reaction rate for 
the same cases. Again, the enhancement of the reaction rate by vortex rollup 
and the influence of the diffusivity are both evident. It is interesting that, 
for the two-dimensional simulations, as the diffusivity is decreased, the peak 
in the integrated reaction rate is delayed somewhat in time. Vortex rollup is 
stretching the reaction interface and bringing the two species in closer 
proximity, but diffusion is still necessary to cause the reaction to occur. 

Also note that, after the rollup has abated (see Figure 23), and the reactants 
in the middle layer have been greatly consumed, the integrated reaction rate 
decreases to near the value for the one-dimensional case. 

The simplest estimate of the average reaction rate RC^C^, the average of 
the product of concentrations, is given by the product of the average concen- 
trations RC^C^. The usefulness of this estimate can be determined by examining 
the inixedness parameter li, defined by 


M - ^5- 

C , c„ 

A B 


(39) 


This quantity is related to the intensity of segregation parameter defined by 
banckwerts (1932) in order to determine the degree of spatial segregation of 
the reacting species. [Related parameters have been defined by, e.g., Sutton 
(1968) and Donaldson and Hilst (1972). j When M approaches 1, the species do 
not tend to be spatially segregated, and the product of the averages RC C is a 
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good estimate of the average reaction rate. However, as M becomes small, the 
reactants are more spatially isolated (at least in the direction of averaging), 
and the product of the averages becomes a poorer estimate of the average 
reaction rate. An integrated mixedness parameter M can be defined by 


M 


C A C B 


(40) 


where the numerator and denominator are appropriate integrated values. 

Figure 32 gives plots of the integrated mixedness parameter versus time for 
the Case 1 simulations and various dif fusivities . Because of the initial con - 
ditions, M is initially 1. As the flow develops in time, M decreases as the 
species become segregated (along the direction of averaging). M also decreases 
with decreasing diffusivity, since sharper reaction zones are defined and the 
species become more segregated. Note that for larger times, after the maximum 
layer extent has been attained and collapse has occurred, M again approaches 1, 
since the reactants become more uniformly mixed in the central region (Figure 25). 
It is clear from these results that the approximation of the average reaction 
rate by the product of the averages would be a poor estimate for these cases, 
overestimating the true average reaction rate by a factor of at least 9. 

It is also of interest to compare the integrated averages for the various 
rollup cases, i.e., for the fundamental mode alone (Case 1), the subharmonic 
mode alone (Case 2), and the fundamental and subharmonic out-of-phase (Case 3). 
Figure 33 contains plots of the integrated products as functions of time for 
these cases. In comparing the results for Cases 1 and 3, we see, somewhat 
surprisingly, that the effect of the second rollup on the integrated product 
is not apparent until after a time of about 22, when the second rollup has 
just about reached its maximum (see Figure 22). Apparently, only when the 
vortex cores begin to merge is the reaction rate enhanced by this second 
rollup. Also, the subharmonic rollup (Case 2) takes more time to enhance the 
integrated product. However, its effect is ultimately greater than the effect 
of the fundamental mode alone, since the subharmonic rollup entrains more 
reactants into the mixing layer. A tentative conclusion from these figures is 
that, for the parameter range considered, the straining of the reaction 
interface is not nearly as important a mechanism in enhancing the reaction 
rates as is the merging of the species in the vortex core. 
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The effect of the different rollup cases on the integrated mixedness 
parameter is shown in Figure 34. We see that the subharmonic rollup and 
double rollup result in even smaller values of M, and that these values are 
maintained for much longer in time, since the collapse and mixing of the 
central layer take much longer to occur. Again, the product of the average 
concentrations would give a poor estimate of the average reaction rate. 

When the integrated product is divided by a characteristic concentration, 
say the free-stream concentration of species A, a length scale results which 
indicates the product thickness. When this length scale is compared to a 
characteristic scale of the velocity field, a nondimens ional estimate of the 
effectiveness of the mixing layer in enhancing the reaction rate is obtained- 
Such a -quantity has been measured experimentally (Breidenthal , 1973; Mungal , 
1983) and is defined by 



where the mean vorticity thickness 6 ^ is defined by 


6 (t) = U 
v 


3z 


U(z , t) 


z=0 


-1 


(42) 


Figure 35 contains plots of this nondimensional product thickness versus time 
for the three different cases computed. Of course, the product thickness is 
initially zero, since the initial product concentration is identically zero. 
The nondimensional thicknesses increase with time until the collapse of the 
layer occurs. After collapse, the product thickness continues to grow, 
whereas the mean vorticity thickness decreases (as does z^, see Figure 22). 
Similar results taken from laboratory 

data and from three-dimensional simulations are discussed in the next section. 

Up to this point we have mainly considered integrated averages, which give 
information about the overall behavior of the reaction. More insight into the 
details of the reaction can be gained by examining spatially averaged 
quantities. For the two-dimensional case, spatial averages are defined [set 
liquation (28)] as an average in the x-direction for fixed z. As was pointed 
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out in Section 4, this spatial average is analogous to a time average in a 
spatially growing mixing layer in the laboratory. 

Figure 36 contains plots of the average concentration of species A, C^, 
versus the transverse coordinate z for several different times from Case 1. 

The growth of the mixing layer in terms of this average concentration is 
apparent. Note that a fairly uniform region forms in the center, which is 
probably an indication of the vortex core. Because of the symmetry in this 
problem, Cg(z) = C^f-z), so that the profile for Cg can be easily obtained from 
that for C^. When this is done, it becomes apparent that there is a large 
region of overlap of and Cg, even though the contour plot (Figure 23) indi- 
cates that the two species are highly segregated. This is the reason for the 
large value of the product of the averages, C A Cg, compared to the average of the 
product, C^Cg. Figure 37 contains plots of the rms values of the concentration 
fluctuations of species A, C^, versus z for the same case and times. Again, 
the growth of the mixing layer is apparent from these plots. When the average 
and rms of the concentrations are compared (Figure 38), it is observed that, 
in the central core and on the upper side of the layer, the rms is almost 
equal to the average. This is a characteristic property of highly intermit- 
tent processes. Finally, Figure 39 contains plots of the average product 
concentration versus z for the same case and times. Note that the product 
field is identically zero at the beginning of the calculation. The growth of 
the product as the reaction proceeds and also the growth of the mixing layer 
are apparent in these plots. The results for the same average quantities 
taken from other two-dimensional cases lend themselves to the same qualitative 
interpretation. 

If we decompose the dependent variables in terms of their average values 
and fluctuations about the average, e.g., 

C = C + C' (43) 

where f' denotes the fluctuation about the average of f, then the average 
reaction rate can be written as 

R Vl ' R I b A f B + (44) 

The term , called the concentration correlation function, is clearly the 

correction needed in the estimate of the average reaction rate by RC.CL. 
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Figure 40 displays these quantities for Case 1 at a time of 12. Note that, as 
indicated by our results for the mixedness parameter, the average product of 
the concentrations is much less than the product of the averages. The correla- 
tion function is almost equal to the negative of C.C , their small difference 
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giving the average reaction rate. Note that many theories suggest that, in 


the fast reaction limit, the sum of C C and C'C' will be of the order of the 

Ad Ad 

second DaiuUohler number taken to the minus 1 power (see, e.g., Gibson and 
Libby, 1972). Also note that a good estimation of the correlation function is 
a necessary ingredient for a useful theory using the moment equations. Again, 
the same statistics were computed for the other two-dimensional cases, ana the 
results were in qualitative agreement with the above. 

Dynamic equations can also be developed for the concentration correlation 
function and for the mean square concentration fluctuations (see, e.g., 

Donaldson and Hilst, 1972). When this is done, the following additional terms 
appear in the equations: C^Cg, C^ 2 , C^ 2 C^, C fi C^Cg , CgC^ 2 , and Cg 2 C^. Note 

that, because of the symmetry of our problem, the first, second, and third 
terms are the symmetric images of the fourth, fifth, and sixth terms, respec- 
tively. In the moment equation approach, the triple correlation terms (as 
well as other turbulent diffusion and concentration dissipation terms) need to 
be modeled. Figure 41 contains plots of C^C^Cg , C^C^, and C^ 2 Cg versus z for 
Case 1 at a time of 12. The term Cg 2 is large and always positive, since 
both and Cg 2 are necessarily nonnegative. On the other hand j- s l ar g e 

and always negative, since is nonnegative and C^Cg k - C^Cg appears to be 
always negative. Finally is roughly asymmetric about the plane of 

symmetry of the mixing layer. The negative values indicate that fluctuations 

in C, are associated with negative values of C., while the opposite holds 

I i A 

true for the positive values- Again, similar results have been obtained for 
the other two-dimensional cases. 

Even more insight into the reaction process can be obtained by examining 
the probability density functions (pdf's) for the species concentrations. 

Because any moment of a random variable can be obtained from its pdf, the pdf 
contains the entire information about these moments. -As discussed in the 
introduction, pdf's have also become an important theoretical/numerical tool 
in the understanding and prediction of reacting flows. Therefore, it is 
useful to examine the behavior of the pdf's for our computed flows. The pdf's 
are computed from our simulations in the following manner. The range over 


7 



which the concentration can vary (in our case from 0 to 1 since the concen- 
trations have been normalized by their free-stream values) is divided into 
equally spaced bins. At a given value of the transverse coordinate z, and for 
every grid point in the x-direction, the value of the concentration is tested, 

and 1/N is added to the appropriate bin, where N is the total number of grid 

points in the x-direction. The resulting value for the n^ 1 bin is just the 
probability for the concentration to be within the limits defined by that bin. 

Figure 42 shows the pdf for the concentration at the transverse 

distance it/ 2 for Case 1 and for various times. We see that for the earliest 

time shown (t = 6), all of the grid points contain concentrations in the range 
of 0.9 to 1, since the reaction had not yet reached this level. The reaction 
zone enters this layer at about a time of 8. At the times of 12 and 18, the 
probability densities have become rather flat, with large values still in the 
bin near 1. This latter fact indicates that an appreciable part of this level 
still contains fluid with only species A. At the time of 24, the probabilities 
of the larger values of concentration have increased considerably, probably 
indicating the collapse of the layer. 

It is also interesting to examine the spatial behavior of the pdf's. To 
do this we choose the time of 12, when the mixing layer for this case has 
reached its maximum extent and has started to collapse. From Figure 23, which 
contains contour plots of for this case and time, we learn that the 
mixing layer has rolled up, and that some fluid with only this species has 
actually penetrated unmixed into the lower region. Strong gradients exist. in 
the braids of the mixing layer, whereas the concentration is more diffuse in 
the vortex core. Figure 43 displays pdf's at four different positions across 
the layer at this same time. At the lowest level (z = — 3 tt / 4 ) , all of the 
concentration values lie near 1, since the reaction zone has not reached this 
level. At the next higher level (z = — TT / 2 ) , the pdf has become flat, with a 
large value in the bin near 1. For the levels nearer the plane of symmetry, 
the profiles are flat, but with large values both near 0 and 1. The large 
value near 0 implies that parts of the fluid at this level contain no species 
A, while the large value near 1 indicates a region of only species A. The 
flat portions of the profiles result from the diffusion and reaction zones. 
These pdf results are entirely consistent with the contour plots and, in some 
ways, help describe them. Again, similar results have been found for other 
two-dimensional cases. 
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6.3 Three-Dimensional Mixing Layer 

The third in our series of simulations was the three-dimensional mixing 
layer. These calculations were intended to simulate three-dimensional 
turbulent flows with chemical reactions, similar to experiments that have been 
carried out in the laboratory. The purposes of the simulations were (i) to 
examine the combined influence of the large-scale structure and smaller-scale 
three-dimensional turbulence on the chemical reactions, (ii) to compare results 
of the simulations with applicable theoretical models, and (iii) to determine 
the validity of the methodology by comparisons with laboratory data. 

The initial conditions for the two concentration fields were the same as 
in the one- and two-dimensional cases [see Equations (17) and (18)]. The 

• ' ' l 

velocity field was initialized using the laboratory data of Wygnanski and 
Fiedler (1970), as discussed in Section 4. A particular realization of the 
initial flow field was intended to model the turbulent flow, and hence was a 
somewhat complex three-dimensional flow. When averages in the x-y plane were 
taken over an initial field [see Equation (28)], some of the resulting average 
properties of the flow were close approximations to corresponding quantities 
taken from the laboratory data, in particular, the mean velocity and turbulence 
intensity profiles and the longitudinal integral scale. Different realiza- 
tions of the statistically same flow field (i.e., taken from the same ensemble) 
could be generated, which give an estimate of the statistical scatter in the 
averaged results. In the following we will often give the results from two 
independent realizations from the same ensemble in order to indicate the 
statistical scatter in the results. 

In the two-dimensional simulations, in order to obtain a single rollup, 
either the fundamental or the subharmonic mode was added to the mean velocity 
to initialize the flow. To obtain multiple rollup, the fundamental and 
subharmonic were added in, but out-of-phase. In the three-dimensional 
simulations, in order to obtain only a single rollup, the extent of the 
computational domain in the x-direction was restricted to be smaller than the 
wavelength of the subharmonic mode (usually this extent was taken to be equal 
to the wavelength of the fundamental mode), so that the subharmonic could not 
develop and a second rollup occur. To compute multiple rollup, the domain 
length was extended to include the subharmonic mode. (Usually, for the 
multiple rollup case, the extent of the computational domain was taken to be 
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the wavelength of the subharmonic mode.) The spectrum of the initial 
turbulence energy [Equation (25)] is very broad, so that appreciable energy 
exists in the fundamental and subharmonic modes to cause the vortex rollup. 

In the three-dimensional simulations, although energy exists at the 
wavelength of the most unstable mode (and at the wavelength of its subharmonic 
when the computational domain is large enough), the most unstable mode itself 
has not been added to the initial field. Also, this field is random and 
somewhat isotropic. Therefore, although vortex rollup is expected to occur, 
it will not take place in as coherent a fashion as in the two-dimensional 
simulations . 

There are a number of ways in which one might generate the initial 
conditions for a turbulent mixing layer, and the choice of the initial 
conditions might have a significant influence in determining the outcome of 
the simulations. For example, it has been found from laboratory experiments 
that the initial conditions can exert a strong influence on the resulting 
mixing layer, especially in the near field of the layer (see, e.g., Hussain 
and Zedan, 1978). Our initial conditions are intended to model a fully 
developed turbulent mixing layer, or the mixing layer produced by fully 
turbulent boundary layers coming off a splitter plate. Such a flow might lack 
the two-dimensional organization which has been observed in many laboratory 
experiments (see, e.g., Browand and Troutt, 1980), although whether the 
large-scale features of turbulent mixing layers are strongly two-dimensional 
is presently very controversial (see, e.g., Dimotakis and Brown, 1976; 
Chandrsuda et al . , 1978). Other possibly reasonable initial conditions 
include (i) the initial mean velocity given by Equation (20) with a low-level 
perturbation field, having a fairly broad spectrum, superimposed, and (ii) 
flows with the same statistical properties as ours in terms of the mean 
velocity, turbulence intensity, and integral scale, but with strong coherence 
in the y-direction, so that the flow will be more well-organized two- 
dimens ionally. Case (i) would involve a calculation of transition to 
turbulence, which is a complete study in itself (see Patera and Orszag, 1981, 
for simulations of a particular aspect of this transition), while the method 
of determining the y-coherence in case (ii) is unclear. Since our results for 
the velocity field were in reasonable agreement with laboratory data using the 
present initialization method (Riley and Metcalfe, 1980a; Metcalfe and Riley, 
1981), we decided to use this method in the present investigation. 
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Figures 44 and 45 give a temporal sequence of contour plots of the 
y-component of the vorticity in two different x-z planes that are located a 
distance /2 apart. We see that initially the vorticity is in a fairly 
narrow region and has a somewhat random orientation. As the flow develops, 
the vorticity begins to become concentrated in local regions, and by a time of 
about 16 there is what appears to be vortex rollup. Note however that rollup 
takes longer to occur than in the two-dimensional simulations (see Figure 19), 
it does not appear to be as strong, the braids between the vortices are not as 
distinct, and there is little lateral coherence in the flow. Figure 46 gives 
a temporal sequence of contour plots of the concentration of species B for 
reactions occurring on the same velocity field and x-z plane as depicted in 
Figure 44. Again, the rollup can be observed. However, the concentration 
gradients do not appear to be as steep as in the two-dimensional case 
(Figure 23), probably due to the small-scale mixing and to the less vigorous 
vortex rollup. 

We next examine the behavior of various length scales that characterize 
the mean velocity field and the mean product field. Figure 47 contains 
results from our simulations for the mean velocity half-width, z M , the 
distance from the plane of symmetry to the transverse plane where the mean 
velocity takes on one-half its free-stream value. Data are presented for two 
different realizations of the single rollup case and tor a double rollup 
case. Similarity theory (see, Tennekes and Lumley, 1972) suggests that this 
quantity, or any characteristic length scale, should vary linearly with time- 
We see that the results are in approximate agreement with this theory. 

Although not shown explicitly in this figure, in the single rollup case, the 
results begin to deviate from linearity at a time of about 24, whereas this 
deviation occurs in the double rollup case for times greater than 36. This is 
expected since, when the second vortex rollup is inhibited by the lack of a 
subharmonic, the layer growth rate is substantially reduced, and in fact it 
can become negative (collapse). These preceding results are consistent with 
those reported by Riley and Metcalfe (1980a; see also Metcalfe and Riley, 
1931). As discussed in Section 3, Riley and Metcalfe (1980a) found that the 
behavior of z,, was consistent with laboratory data if space and time were 
assumed to be related by the transformation given by 

x = | (U L + U 2 )t 


41 



The present results are also consistent with laboratory data if the same 
assumption is made. Figure 48 contains similar data for the mean vorticity 
thickness [see Equation (42)]. We see that again the growth rate is 
approximately linear, consistent with similarity theory. 

The mean product thickness, Zp, defined by 




C p (£,t) d£ 


( 45 ) 


was also computed from our simulations. Figure 49 contains plots of the mean 
product thickness versus time. Included are data from two different realiza- 
tions for the single rollup case, from a double rollup case, from a simulation 
with the molecular diffusivity D increased by a factor of 2 over its usual 
value, and for an infinite reaction case. Again, similarity theory suggests 
that this length scale should increase linearly with time, and again the 
results are in approximate agreement with the theory. Comparing the results 
for the two independent realizations, we find that the results deviate by 
about 5% to 10% at the end of the calculation. A comparison of the results 
for the two different dif fusivities indicates that increasing the diffusivity 
by a factor of 2 results in a less than 10% increase in the product thickness. 
In contrast, the theoretical result of Burke and Schumann (1928) for the 
interdiffusion and reaction of two species with no velocity field predicts an 
increase of about 40%. Comparing the single and double rollup cases leads to 
the conclusion that the results are approximately the same up to a time of 
about 20, when the growth in the single rollup case is inhibited by the lack 
of a subharmonic. Finally, comparing the results for the infinite reaction 
case with those for the finite reaction cases, we find that the infinite 
reaction rate leads to about a 15% increase in the product thickness. 

As discussed in the previous subsection, the mean product thickness 
divided by the mean vorticity thickness is a useful parameter with which to 
compare experimental data. Figure 50 is a plot, taken from Mungal (1983), of 
this thickness ratio as a function of the equivalence ratio (<j> = C^/Cg^). 
These data were taken from experiments performed in a gaseous mixing layer, so 
the Schmidt number was approximately 1 as in our simulations. For an 
equivalence ratio equal to 1, which is then the same as in our simulations, 
Mungal's thickness ratio equals approximately 0.35. Similar data from Wallace 
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(1981) give a value of about 0.3. Breidenthal (1978) measured the thickness 
ratio in experiments in water. His Schmidt number was approximately 600, and 
his equivalence ratio was about 200. If the experiments in water have 
approximately the same dependence ou equivalence ratio as the previously 
discussed experiments in a gas, then Breidenthal would obtain a thickness 
ratio of about 0.23 for an equivalence ratio of 1. A recent model by 
Broadwell and Breidenthal (1982) suggests that the differences between the 
results of Breidenthal and those of Mungal and Wallace are due to the 
substantial difference in the Schmidt number. 

Figure 51 contains results for the thickness ratio as a function of time 
from our various three-dimensional simulations. We see that the thickness 
ratio starts from zero and ultimately attains an approximately constant value 
of about 0.2 for the finite reaction rate case and 0.22 for the infinite 
reaction rate case. Thus, the results are very close to those of Breidenthal, 
but are about 25% to 30% below those of Mungal and Wallace. Since the value 
of the Schmidt number in our calculations is 1, we would expect better 
agreement with the latter experiments. However, the above agreement must be 
considered reasonably good, especially since there are no adjustable 
parameters in our calculations. 

Another nondimensiona! quantity of interest is the integrated mixedness, 
defined by Equation (40). Figure 52 contains results for this quantity 
computed from various three-dimensional simulations. Comparing these results 
with similar results from two-dimensional simulations (Figure 32), we see that 
the mixedness in the former case is not as small as in the latter, indicating 
that the reacting species are not: quite as segregated in the three-dimensional 
case as in the two-dimensional case. This is probably due to the small-scale 
three-dimensional turbulence and also the lack of coherent vortex rollup in 
the three-dimensional simulations. 

We next examine the behavior of some of the spatially averaged quantities 
[see Equation (28) for a definition of the spatial average]. Figure 53 contains 
a sequence of plots of the mean concentration of species A, C^, for a number of 
selected times for the first realization of the double rollup case. The growth 
of the mixing layer and the depletion of this species by the reaction are 
evident from this plot. Figure 54 gives a similar sequence for the concentra- 
tion of species B for the same conditions. Note that, due to the symmetry of 
this problem, Cg(z) = should hold. This is approximately the case, and 
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the deviation from it gives some indication of the statistical scatter in our 
results. Note also the large overlap in the profiles for and Cg for the 

later times. Figure 55 contains a sequence of plots for taken from the 

second realization at the same times. As described at the beginning of this 
section, comparing the results again gives a good indication of the 
statistical scatter in our results. 

Figure 56 contains plots of the rms of the fluctuating concentration of 
species A for the first realization and the same times. Initially, this 
fluctuating concentration is zero. However, it rapidly grows and attains an 
approximately constant maximum value of about 0.31 C . The profiles tend to 

Aco 

be somewhat skewed to the lower side of the mixing layer, since this is the 
region containing species A. The growth of the layer is also apparent from 
this plot. Figure 57 provides the same profiles for species B for the same 
conditions. Again the symmetry condition C_(z) = C.(-z) should hold, which is 
approximately the case. Figure 58 contains similar profiles for species A 
taken from the second realization. Again, there is approximate agreement 
between the results for the two realizations. 

Figure 59 contains plots of RC.C , RC.Ci, and RC'C' for the first 
° r A B’ A B A B 

realization at the time 12. From Equation (44) we know that the first term, 
the mean reaction rate, is equal to the sum of the second two terms. It is 
clear from the figure that the product of the means is, as in the 
two-dimensional case, not a good estimate of the mean of the products. These 


results are typical of our other three-dimensional calculations. 

Figure 60 contains plots of C. C'Ci, C. C ' 2 , and cT 2 cT for the first 

A A d A d Ad 

realization of the time of 12. From symmetry conditions, these should be 
symmetrically related to Cg C^Cg, Cg C^, and Cg Z C^, respectively. Note that, 
as discussed in the previous subsection, these terms are of importance in the 


equations for C' 2 , C ' 2 , and C'Cl 
_ A B A_B 

C A is large and positive, 


As in the two-dimensional case, the function 
C^Cg is large and negative, and Cg^C^ is 


roughly asymmetric about the axis of symmetry. The interpretation of these 


results is qualitatively the same as for the two-dimensional case. Apparently, 


the two-dimensional dynamics govern the basic behavior of these quantities. 


Figure 61 contains plots of the same quantities for the second realization at 
the same time. Note the significant scatter in the result for C A 2 Cg. It is 
common in analyzing statistical data that the errors increase with the order 


(power) of the statistical quantity. 
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A theoretically and experimentally useful method to analyze results for 
flows with simple geometries is in terms of similarity theory. Similarity 
theory implies that there is a characteristic length scale (&), velocity 
scale (U, and hence time scale H/U), and concentration scale (C) which 
should characterize all of the results of an experiment. (See, e.g., Tennekes 
and Lumley, 1972, for a lengthy discussion of similarity theory.) For a 
temporally developing flow, these scales should only depend on time, and, for 
example, the mean velocity should be given by 

u(z,t)/U(t) = f(z,t) = g[z/£(t)] (46) 

Therefore, when similarity holds, if the mean velocity is normalized by U and 
the transverse coordinate z by i, then the data should collapse when plotted 
for various times. In the case of the mixing layer, the characteristic 
velocity scale should be a constant in time and can be taken to be the 
velocity difference across the layer, U. The characteristic concentration 
scale should also be a constant, and can be taken to be the free-stream value 
of one of the reacting species, e.g., C Aoq . The length scale should vary 
linearly with time and can be taken to be the mean velocity half-width or 
the mean vorticity thickness 6 y , which should be proportional. It has 
been found from laboratory experiments for turbulent mixing layers that the 
statistics of both the velocity field (see, e.g., Wygnanski and Fiedler, 1970; 
Brown and Roshko, 1974) and the reacting species fields (see, e.g., Konrad, 

1977; Mungal , 1983; Batt, 1977) agree well with similarity theory. Because of 
the theoretical implications of this theory, and also because of the good 
agreement of this theory with laboratory data, it is important to determine 
the consistency of our computed results with this theory. 

We have already presented data for the time dependence of various length 
scales in our computed flows. We found that both the mean velocity half-width 
(Figure 47) and the mean vorticity thickness (Figure 48) vary approximately 
linearly with time, in agreement with similarity theory. We also found that 
the product thickness (Figure 49) grew approximately linearly with time, again 
consistent with the theory. Thus, all three scales remain proportional, and a 
single scale characterizes the flow field. We have found from previous 
studies (Riley and Metcalfe, 1980a; Metcalfe and Riley, 1981) that the mean 
velocity, turbulence intensity, and Reynolds stress were all in approximate 
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agreement with similarity theory. We next explore whether the results for the 
concentration fields are consistent with this theory. 

Figure 62 contains a sequence of plots of the average concentration of 
species A plotted in similarity coordinates for various times in a double 
rollup, finite reaction rate calculation. The concentrations have been scaled 
by the free-stream value of C^, while the transverse distance z has been 
scaled by the mean velocity half-width z^, computed from the mean velocity 
profiles. We see that, except for the initial condition, the collapse of the 
data is excellent. Note that the width of the mixing layer as measured by 
z^ grows by a factor of more than 5 over the course of this calculation, so 
that the average concentration itself has changed significantly over this time 
frame (see Figure 53). Figure 63 is a similar plot for the concentration of 
species B. Again the collapse of the data using the similarity scaling is 
very good. Also note that the symmetry in the problem implies that 
Cg(z,t) = C A (-z,t). Comparing Figures 62 and 63, we see that this is approxi- 
mately the case. 

Figure 64 contains a similar sequence of plots of the rms concentration 
fluctuation of species A taken from the double rollup results. Initially, the 
concentration fluctuations are zero, due to the initial conditions. However, 
the concentration fluctuations adjust and become approximately self-similar 
with a peak value of about 0.29. The profiles tend to be somewhat skewed to 
the bottom of the mixing layer, the region containing species A. This 
behavior is qualitatively similar to that for laboratory data (see, e.g., 
Konrad, 1977). Some insight into these results can be obtained by examining 
the dynamic equation for the mean square concentration fluctuation, i.e., 
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(47) 


and noting the behavior of the various concentration correlations given above 
(see Figure 60). We find that the principal generation terms appear to be 
W '^A l$z ^A’ can interpreted as the generation of concentration 

fluctuations resulting from turbulent diffusion along the mean concentration 
gradient (a work-like term if the analogy between concentration fluctuations 
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and turbulence kinetic energy is made), and RC^ C^Cg, which also appears as a 
source term in the equation for C^ 2 . The term w' C^ 2 / 2 represents the turbu- 
lent diffusion of C ' 2 . C ' 2 is dissipated by the last term on the right hand 

A A 

side of Equation (47), the usual turbulence dissipation term, and RCj^2 Cg, which 
represents the loss of C^2 resulting from the decrease in due to the chemical 
reaction. Not that the peak in C'2 tends to lie near the high gradient region 

A 

of (compare Figures 62 and 64). This indicates that the work-like term, 
which generally has its maximum near the region with the largest gradient in C^, 
is a significant contributor to the generation of C^ z . We hope to further 
investigate the dynamics of the concentration fluctuations in the future by 
carefully comparing the various terms in Equation (47). Figure 65 contains 
similar plots for the concentration fluctuation for species B. Again, the 
agreement with similarity theory is rather good. From the problem of symmetry, 
the relation C (z,t) = C (-z,t) should hold, and comparing Figures 64 and 65 we 

D A 

find that this is approximately the case. 

We next examine the behavior of the concentration fluctuation correlation, 
which is plotted in Figure 66 in similarity coordinates for a sequence of 
times for the double rollup case. This quantity is also initially zero 
because of the initial conditions. However, it also becomes approximately 
self-similar, with a peak value of about -0.03. There is somewhat more 
scatter in this result, which is expected since the statistical error in 
computing correlation functions is generally higher than when computing 
averages and mean squares. Symmetry conditions imply that this correlation 
function should be an even function of z, which is approximately the case. 

It is interesting to examine the effect of the reaction rate on the 
concentration statistics. Figure 67 contains similarity plots of the average 
concentration of species A for various times taken from a double rollup, 
infinite reaction rate calculation. The velocity field used in these 
calculations was identical to that used in the previous finite rate 
calculations. We see that the similarity scaling collapses the results very 
well for this case, and, comparing Figures 62 and 67, that the reaction rate 
has very little effect on the average concentration. Figure 68 contains a 
similar plot for the rms concentration fluctuation for the infinite reaction 
rate case. Again, agreement with similarity scaling is rather good. Comparing 
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Figures 64 and 68 we see that the rms concentration is about 15% higher in the 
infinite reaction rate case. 

For the infinite reaction rate case, the product field was also computed. 
Figure 69 contains a sequence of similarity plots of the average concentration 
of the product for various times taken from the same double rollup calculation 
as reported above. Initially, the product field is finite and is determined 
by both the initial conditions on the concentration fields [Equations (17) and 
(18)] and also the interpretation of the product field in terms of conserved 
scalars [see Equation (16) and following]. The product field initially 
adjusts and then shows good agreement with similarity scaling, attaining a 
peak value of approximately 0.30. The average product profile has an 
approximately symmetric, Gaussian shape, which is qualitatively similar to 
laboratory results for experiments with equivalence ratios of about 1 (see 
Mungal , 1983). 

Figure 70 contains similar plots for the product concentration fluctuation. 
This quantity is initially zero and also becomes approximately self-similar, 
with peak values of roughly 0.15. Note the double-peaked behavior of the 
plots, which is consistent with laboratory data (Konrad, 1977; Mungal, 1983). 
This behavior is probably due to the generation of concentration fluctuations 
occurring in the region of large gradients of Cp, and is reflected in the 
turbulence work-like termw'C^ Cp in the equation for C^ 2 . 

We find that the results of our simulations allow extensive examination 
and interpretation, both in terms of contour plots of relevant variables and 
also in terms of various statistical quantities of interest. The results also 
show good agreement with similarity theory, as does the laboratory data. This 
is true even for quantities that are initially far from their similarity 
values. And, without any adjustment of parameters, the results for the 
nondimens ional product thickness are in rough agreement with laboratory data. 
We next briefly examine some of the implications on existing theories of the 
results of our calculations. 

One modeling approach using the Reynolds-averaged equations (Donaldson and 
Hilst, 1972) includes, in addition to the equations for and Cg, also the 
equations for C^ 2 , Cg 2 , and C^Cg. These equations can be closed by assuming 
that triple moments such as Cj^Cg are small compared to lower order terms like 
C A C A^B anc * C A 2 t see E 9 u ation (47)]. However our results indicate that the 
triple moments are not small in comparison to the lower order terms (see 
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Figure 60), indicating that a model based upon this assumption will probably 
lead to poor, predictions.. 

Another approach to closing the Reynolds-averaged equations is due to 
Mason and Spaulding (1973), and has been termed the eddy-breakup model. A 
version of this model which is currently being used by Sturgess (1983) assumes 
that the average reaction rate term is given by 

RC.C = min(C C )/T (48) 

A D A D 


f C. and 
A 

reactant, and T is a turbulence time scale. Similarity theory predicts that T 
should be constant across the mixing layer and should increase linearly with 
time. Figure 71 is a plot of the mean reaction rate divided by the eddy- 
breakup model prediction as a function of the transverse coordinate z. [No 
attempt has been made to set the proportionality constant implied by 
Equation (48).] This result is taken from a double rollup, finite reaction 
rate simulation at a time of 12 but is typical of results computed at other 
times and from other simulations. The breakup model predicts that this ratio 
is constant across the layer. However, the value varies appreciably across the 
layer, indicating a significant difference between the theory and our calcula- 
tions. Figure 72 displays a plot of the ratio min(C A> C^/C^C^, computed along 
the centerline of the mixing layer, versus time. Again, the result is taken 
from the double rollup, finite reaction rate case but is consistent with other 
simulations. The breakup model predicts that this ratio should vary linearly 
with time, which is in approximate agreement with our results. 

A critical term in closing the equations for the average concentrations is 
the correlation For statistically homogeneous flows, Toor (1969) made 

assumptions about the probability distributions of the reacting species and 
found that this correlation was identical for the very fast and very slow 
reaction limits. He then conjectured that the correlation would be the same 
for all reaction rates. If this is truly the case, then the correlation for 
the reacting flow case could be predicted using the correlation for the non- 
reacting flow case. The latter is much easier to predict. It is of interest 
to see how well this hypothesis holds up for our (nonhomogeneous ) simulations. 
In addition to simulations for the finite and infinite reaction rate cases, we 
have also performed simulations for the limiting case of no reaction. 


CL, i.e., the concentration of the lean 


Here min(C, 
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Figure 73 contains plots of the concentration fluctuation correlation 
taken from two-dimensional, single rollup simulations for a nonreacting case 
and for a finite reaction case. We see that the agreement between the 
correlations taken from these two cases is quite good. This result is typical 
of comparisons taken from the two-dimensional simulations. Figure 74 shows a 
similar comparison taken from three-dimensional simulations. We see that the 
agreement is not as good as for the two-dimensional case. However, the 
results do suggest that the theory still might be useful in predicting the 
correlation. 
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7. CONCLUSIONS AND DISCUSSION 

The work discussed in this report has involved the following tasks. First 
of all, two computer codes were developed, a two-dimension plus time computer 
code and a three-dimension plus time computer code. These codes solve the 
dynamic equations for th£ diffusion-reaction problem [Equations (l) to (5)] in 
two and three dimensions, respectively. Next, these codes were extensively 
tested to determine the reliability of the numerical methods and the parameter 
ranges over which accurate solutions could be expected. Finally, the 
following sequence of problems was computed: (i) reactions on a unidirection; ; 

(one- dimensional) mixing layer, (ii) reactions on a mixing layer experiencing 
large-scale, two-dimensional vortex rollup, and (iii) reactions on a three- 
dimensional turbulent mixing layer. 

The numerical testing involved the comparison of computed results with ! 

i 

exact solutions for a number of different cases. Both rigid body rotation and 
vortex rollup flow fields were used. The work greatly extended the results of 
Oirszag (1971) for the advection of a passive scalar on a rigidly rotating floto 
field (the color problem) to include also diffusion, chemical reaction, and 
more complex flows. We have found that high accuracy of the spectral methods 

; i 

observed in the advection case is also obtained when these further compli- 
cations are present. For computing reacting species, it was found useful to 
introduce a weak spatial filter [Equation (11)] to insure numerical 
stability. Our results indicate that spectral numerical methods may prove to 
be useful in the future both for solving the model equations for combustor 
processes as well as for future studies of chemically reacting turbulent flows 
employing direct numerical simulations. 

The approach of direct numerical simulations allowed extensive examination 
and interpretation of the reaction process. From the one-dimensional 
simulations, we found that we could easily compute finite reaction rates near 
the fast reaction limit, that the results were fairly insensitive to the 
initial conditions (for the class of initial conditions computed), and that 
the results were in reasonable agreement with theoretical predictions. 

For the two-dimensional simulations, we first examined contour plots of 
the vorticity and concentration fields. We found that vortex rollup both 
caused large regions of fluid containing one species to be drawn deeply into 
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the regions originally containing the other fluid and also significantly 
stretched the reaction zone. However flame shortening, which to some extent 
moderated this latter effect, was noticed in the cores of the vortices. In 
all of the cases computed, the vorticity field and the product field 
approximately coincided. 

From the computation of volume averages, we observed the enhancement, o f 
the overall reaction rate due to the vortex rollup. It also appeared that the 
merging of vortex cores was a more significant mechanism in increasing the 
overall reaction rate than the straining of the reaction interface. And the 
value of the nondimens ional product thickness, the ratio of the product 
thickness to the mean vorticity thickness, was found to be about 0.07., From 
the spatially averaged (in the x-direction) data, the growth of the mixing 
layer, the depletion of the reactant, and the creation of products were 
observed. Significant species segregation was apparent, so that, for example, 
the average of the product of the reactant concentrations was very different 
from the product of the averages. Instead, the product of the average 
concentrations was approximately equal to and opposite the correlation between 
the fluctuating concentrations. Some of the higher order correlations were 
also examined. Probability density functions were computed, which were found 
to be useful in interpreting the various concentration fields. 

In the three-dimensional simulations, the contour plots indicated that the 
vortex rollup takes longer to develop, and the vortices and braids are not as 
distinct as in the two-dimensional case. Also, the vortices that develop are 
not strongly correlated laterally. Both the contour plots and the statistical 
results indicated that the spatial segregation was also not as strong as in 
the two-dimensional case, probably due to the weaker vortex rollup as well as 
the effects of smaller-scale , three-dimensional turbulence. 

Comparisons were made between the simulation results and results using 
similarity theory. Approximately linear growth rates of various computed 
length scales, including the mean velocity half-width, the mean vorticity 
thickness, and the mean product thickness, were obtained and were in agreement 
with the theory. Similarity scaling was found to collapse quite well the 
results for the average reactant concentrations, the rms fluctuating reactant 
concentrations, the concentration correlations, the average product 
concentrations, and the rms fluctuating product concentrations. This collapse 
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was attained over a time interval in which the width of the layer grew by a 
factor of over 5, even though the initial conditions for these quantities 
differed significantly from the similarity values. 

Some limited comparisons were made with laboratory data. Of course, the 
agreement with similarity theory is consistent with the laboratory results, 
since the laboratory data also agree very well with the theory. Computed 
profiles that were qualitatively similar to corresponding laboratory profiles 
were obtained for the average reactant concentrations, the rms fluctuating 
reactant concentrations, the average product concentrations, the rms 
fluctuating product concentrations, and the concentration correlations. 
Finally, the nondimens ional thickness ratio was computed to be approximately 
0.22 for the infinite reaction case. This is to be compared to the value of 
0.3 and 0.35 computed from the data of Wallace (1981) and Mungal (1983), 
respectively, for flows with Schmidt numbers of order 1, and 0.23 estimated 
from the data of Breidenthal (1978) for a flow with a very large Schmidt 
number. Since our Schmidt number is equal to 1, we would expect better 
agreement with the former results. The agreement must be considered 
reasonably good, since there were no adjustable parameters in the simulations. 

We have made some comparisons with existing theories. Donaldson and Hilst 
(1972) have suggested, in addition to using the equations for the average 
concentrations, including the equations for the concentration fluctuations and 
correlation. These equations can be closed by neglecting certain triple 
moments when compared to certain lower order terms. However, our results 
indicate that the triple moment terms are as important as other terms in the 
equations, so that an assumption of this type will probably lead to poor 
predictions. Mason and Spaulding (1973) have proposed a model for the mean 
reaction term, suggesting that it will be proportional to the average 
concentration of the lean species divided by a turbulent time scale. We have 
found that such an assumption will only be moderately successful if applied to 
our case. Finally Toor (1969) has suggested estimating the concentration 
correlation of the reacting species in terms of that for the nonreacting case, 
which is much easier to model. Although this was proposed mainly for 
statistically homogeneous flows, we find that it is a reasonable approximation 
for our reacting flow simulations. 
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While our present results inspire optimism in the application of the 
techniques of direct numerical simulations to chemically reacting turbulent 
flows, one should still keep in mind the potential differences between our 
computed flow fields and laboratory experiments. First of all, our simulations 
have been carried out at a rather low turbulent Reynolds number [R^ of 
about 50 compared with values of several hundred for the laboratory experi- 
ments (see, e.g., Mungal , 1983)]. Since the Reynolds number is significantly 
lower, the rate of fluid surface (and hence reaction interface) stretching, 
which depends directly on the small-scale vorticity field, is diminished. 

This latter fact implies that the reaction interface does not increase as 
rapidly in our simulations as in the laboratory experiments, thus possibly 
decreasing the overall reaction rate. Most theories imply, however, that if 
the Reynolds number is large enough, the overall reaction rate should not 
depend on the Reynolds number but only on the rate at which the large-scale 
structures bring the two species together. The recent theory by Broadwell and 
Breidenthal (1982) suggests in addition that the overall reaction rate should 
also depend on the Schmidt number. The question remains as to whether the 
Reynolds number in our simulations is high enough for Reynolds number 
similarity to be valid. Previous computations of (nonreacting) wakes (Riley 
and Metcalfe, 1980b) and mixing layers (Riley and Metcalfe, 1980a; Metcalfe 
and Riley, 1981) indicated that Reynolds number similarity is at least 
approximately attained in these simulations. However, a higher Reynolds 
number may be needed for similarity in reacting flows. 

It must also be remembered that our calculations are for the temporally 
growing mixing layer, whereas all of the reacting flow experiments have been 
for a spatially growing layer. This difference can be exemplified by the 
entrainment ratio, defined by Konrad (1977) to be the average amount of fluid 
entrained into the mixing layer from the high-speed side divided by the 
average amount of fluid entrained from the low-speed side. For constant 
density flow with a free-stream velocity ratio of 0.38, Konrad estimated that 
the entrainment ratio was approximately 4:3, so that 4 parts of fluid from the 
high-speed layer are entrained into the mixing layer for every 3 parts of 
fluid from the low-speed layer. Of course, in the temporally growing layer, 
because of the symmetry in the problem, the entrainment ratio is 1. 
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The initial (upstream) conditions for our flow fields may be considerably 
different from those in the reported laboratory experiments, perhaps causing a 
significantly different flow field to develop. Our initial flow is rather 
randomly distributed in space (see Figures 44 and 45), whereas in the 
laboratory, the flow just downstream of the splitter plate is somewhat 
two-dimensional, even if the boundary layer on the plate is very turbulent, 
since the splitter plate itself is two-dimensional. It has been established 
that initial conditions in the mixing layer can have a significant effect on 
the subsequent layer development, especially in the near field (see, e.g., 
Hussain and Zedan, 1978). This effect can become more pronounced as more 
sensitive statistical quantities are examined. The comparisons of our computed 
results with laboratory data for the case of the mean velocity, turbulence 
intensities, and Reynolds stress (see Riley and Metcalfe, 1980a; Metcalfe and 
Riley, 1981) have indicated that the initial conditions were adequate to model 
the essential physics leading to these quantities. However, it is possible 
that the reaction rates are more sensitive to the initial conditions. 

In order to further explore the effects of initial conditions, simulations 
with different initial conditions should be computed. Two possible candidates, 
as discussed in Section 6.3, are the following: (i) the initial mean velocity 

given by Equation (20) with superposition of a low-level perturbation field, 
which has a relatively broad spectrum, and (ii) flows with the same 
statistical properties as ours, but with stronger coherence in the lateral 
direction, so that the flows will be more well-organized in that direction. 
Based upon laboratory experiments (Konrad, 1977), case (i) should produce 
vortex rollup followed by a three-dimensional breakdown consisting of 
longitudinally oriented vortices. Some simulations of this type (for 
nonreacting flows) have been carried out by Patera and Orszag (1981). For 
case (ii), the lateral coherence introduced should be consistent with the 
results on lateral coherence presented by Browand and Troutt (1980). 

In addition to further work on initial conditions, more comparisons with 
laboratory data are needed in order to determine the level of confidence of 
the simulation methodology. In particular, detailed comparisons of mean and 
fluctuating concentration profiles (for both the reactant and the products) 
should be made. During the course of this study, data could not be found with 
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which to make such comparisons. However, it appears that if the data of 
Konrad (1977) and Mungal (1983) are analyzed properly, then such comparisons 
can be made. Also, some ongoing experiments by Choudbury et al . (1983) may 
produce results against which to compare our numerical simulations. 

Additional simulations can also be performed to learn more about the 
effects of turbulence on the reaction process. For example, the dynamics of 
the reaction can be examined in more detail by computing the various terms in 
the appropriate dynamic equations [for example in Equation (47)] and 
determining which of the terms are dominant and, thus, which effects are of 
importance and should be included in any model. Furthermore, calculations can 
be performed for different choices of parameters. For example, the effects of 
differences in the diffusivities of the two reacting species could be 
explored, problems with nonstoichiometr ic free-stream conditions could be 
studied, and cases with low- to moderate-speed reactions could be computed. 

Of additional importance would be to examine cases with more complex 
reactions, possibly with several reactants and the reactions proceeding in 
both directions. Use could be made of the scalar approach to compute cases 
where equilibrium chemistry is assumed (Bilger, 1980). 

Finally, of utmost importance is to relax the assumption of no heat 
release and compute cases where the heat release is allowed to effect both the 
flow field and the reaction rates. This should probably initially be done in 
the low Mach number limit so that acoustic waves can be filtered out of the 
problem (see, e.g., Oran and Boris, 1981). The advantages of this approach 
are that the numerical time stepping is then not controlled by the acoustic 
waves, so that the time steps can be increased significantly, and that no 
radiation boundary conditions are needed. However, the basic physics of 
interest, the effects of heat release on the turbulence and reaction rates, is 
retained . 
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APPENDIX A 


NUMERICAL RESOLUTION PROBLEMS 


In reacting flows of the type addressed in this report, i.e., mixing 
layers, that involve the interdiffusion and reaction of nonpremixed species 
and in which the reaction rate is fast relative to the species diffusion, i.e., 


2 

RL Z C 


>> 


1 


(see Table 1 for definition of the parameters), very steep gradients can 
develop in the concentration fields, even from smooth initial conditions. Of 
course, in the limit of an infinite reaction rate, the reacting species become 
segregated, and the spatial derivatives of the concentrations are discontinuous 
at the interface between the species. These steep gradients cause several 
numerical problems. First of all, the number of modes or grid points 
necessary to resolve a concentration field to a certain level of accuracy 
increases with the steepness of the gradients in the field (assuming other 
factors, such as the size of the computational domain, remain constant). 

Thus, a calculation can be started with adequate resolution for all the 
initial fields, but can quickly break down as steep gradients develop. Since 
the extent to which this happens can depend strongly on the evolution of the 
flow field, its occurrence can be very difficult to predict a priori- 
Therefore, a number of calculations may have to be made with increasing 
spatial resolution to determine if an accurate calculation is possible for a 
particular flow and species field configuration with a given choice of 
parameters (e.g., R and D). 

An attempt to resolve concentration fields with steep gradients can lead 
to a second problem, which can also have a detrimental effect on the accuracy 
of the calculation. Near a region where the concentration field increases 
rapidly from a value of zero, the numerical approximation can produce a 
"ringing" or Gibbs phenomenon effect (cf. Wengle and Seinfeld, 1978; McRae 
et al . , 1982) that results in spurious negative values of C^. As can be 
seen by an examination of the diffusion-reaction equations [Equations (l) and 
(2)], if this occurs where is positive, as it might near a reaction zone, 
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the net effect is an artificial production of C g and an artificial destruc- 
tion of the product. While this effect may be somewhat mitigated by the 
tendancy of the C A equation to drive C A back toward zero, if the steep 
gradient and hence the negative overshoot of C A persists, a significant 
error can accumulate over a long time integration. The problem becomes more 
serious if spurious negative regions of C and C D overlap. For example, 
neglecting advection and diffusion, and assuming C A = C g < 0 and R = 1, 
Equation (1) behaves like 



C A " (C A t + 1) 
o 

where C is the initial value of C A . This becomes singular when 


(Al) 


(A2) 


t 



o 

Several approaches have been suggested to address this problem. Cullen (1976) 
and Raymond and Gardner (1976) have added diffusive or dissipative terms to 
the equations, and Storch (1978) has added high wave number filtering to damp 
out the oscillations. Forester (1977) has devised a spatially local nonlinear 
filtering scheme to decrease the amplitude of the oscillations. A number of 
other alternatives are discussed by McRae et al. (1982). The most direct 
approach is to simply enforce the condition [Equation (11)] 


C (x , t ) 


C(x , t ) 
0 


C > 0 
C < 0 


(A3) 


where C is the value of the concentration prior to applying the restriction. 
Although this has the apparent disadvantage of violating mass conservation, to 
the extent that the calculation is accurately resolving the steep concentration 
gradients, the net effect should be small. In fact, if the results of the 
simulation are significantly dependent on the filtering or smoothing scheme 
chosen, the spatial resolution in the calculation may well be inadequate. 
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In our simulations, we have found that the application of condition (A3) 
every 5 to 10 time steps works very well in preventing the solutions from 
diverging, and that the results are quite insensitive to whether the condition 
is applied every time step or every 10 steps. This suggests that, with the 
high accuracy of the spectral methods we have employed and in the parameter 
ranges we have chosen, the approximation errors introduced by the Gibbs 
phenomenon are small. Also, the most important objective under these 
conditions is to employ a simple scheme that prevents these small errors from 
amplifying over the course of the simulation. We have found that our 3D 
simulations with a turbulent velocity field are less sensitive to the effects 
of oscillation errors near the interface than are the 2D simulations. For 
example, we have performed run JN20B both with and without the application of 
condition (A3) and the energies and other moments differ by less than 0.1% 
between the two runs. A possible explanation for this is that in a highly 
structured, coherent flow field like that in our 2D mixing layer simulations, 
it is possible for regions of overlapping spurious negative concentrations to 
persist for some time and thus cause the calculation to diverge. In the 3D 
simulations, however, such regions and the steep gradients that generate them 
are being advected and deformed much more drastically by the turbulent flow 
field. The turbulent diffusion acts both to decrease the steepness of the 
concentration gradients and to prevent possible negative concentration regions 
from overlapping for a sufficiently long time interval to destabilize the 
calculation . 


59 


APPENDIX B 


THE RELATIONSHIP BETWEEN MONTE-CARLO METHODS 
AND DIRECT NUMERICAL SIMULATIONS 


The basic ideas- for the application of Monte-Carlo methods to probability 
density function (pdf) equations stem from the theory of-Markov processes and, 
in particular, Brownian motion (see, e.g., Wax, 1954). A dynamic system is 
assumed to be described by the following equation: 


dv 

dt 


+ v 


f ' 


(Bl) 


where v(t) is, for example, the speed of a particle, and f(t) is a forqing 
function. Here, f is assumed to be a known white noise process, so that v is • 
a stochastic process. The equation for the pdf of the displacement x can be 
found to be : 


3P _ 9jP 

3t " 3x 2 


(B2) 


This equation is often referred to as the Fokker-Planck equation, and the 
dynamic equation is called the Langevin equation. If one is interested in 
solving an equation of the form given by Equation (B2), then the Monte-Carlo 
approach is to numerically solve Equation (Bl) a large number of times for 
random, independently selected initial conditions, and to compute the pdf or 
its moments by averaging over the ensemble of solutions. Note that Equation 
(Bl) is of lower dimension than Equation (B2) and, thus, is generally 
significantly easier to solve numerically. 

Applying these concepts to the pdf equations for reacting turbulent flows, 
the pdf equation is analogous to the Fokker-Planck equation. What is then 
needed is an analogous Langevin system. In coalescence-dispersion modeling, 
this system was invented by Curl (1963), and Monte-Carlo methods have been 
used to solve for the pdfs and their moments (e.g., Spielman and Levenspiel, 
1965; Pratt, 1976). Several authors (see, e.g., the review by O'Brien, 1981) 
have derived more rigorous pdf equations describing reacting turbulent flows, 
although adhoc assumptions are s til 1, required . Pope (1979) has devised an 
analogous system for one of these approximate pdf equations, making the 
application of Monte-Carlo methods again possible. 
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Of course, the analogous Langevin system for the exact (unmodeled) pdf 
equations is just the coupled Nav ier-Stokes and diffusion-reaction equations. 

A Monte-Carlo solution to the exact pdf equations would require solving this 
coupled system of equations many times and taking averages over the computed 
ensemble. However, for a dynamic system that is statistically homogeneous in 
either space or time, spatial or temporal averages are, under usual condi- 
tions, equivalent to ensemble averages. Thus, direct numerical simulations 
can be considered as a type of Monte-Carlo method, where spatial (or temporal) 
averages, sometimes supplemented by ensemble averaging over a small number of 
realizations, are used instead of ensemble averages. Direct numerical 
simulations have the advantage of a more exact Langevin system than the 
approaches of Curl or Pope, so that the physical system is modeled better. 
However, they have the disadvantage of a much more difficult numerical 
calculation than the previous approaches, so that, in particular, they are 
more difficult to apply to complex problems. 
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LIST OF SYMBOLS 

molar concentration of species i 
free-stream value of C. 

l 

molar concentration of product 
reactant diffusivity 
first Damkohler number 
second Damkohler number 
species conservation invariant 
truncation wave number vector 
wave number vector 

length nondimensionalization factor = UT 

length of computational domain in streamwise (x) direction 
length of computational domain in lateral (y) direction 
length of computational domain in transverse (z ) direction 
mixedness parameter 

number of Fourier modes in each spatial direction 

stoichiometric coefficient 

number of grid points in x-direction 

number of grid points in y-direction 

pressure 

reaction rate constant 
instantaneous reaction rate 
Reynolds number 
Reynolds number = UL/v 
turbulent Reynolds number = uA/v 
Schmidt number 

time nondimensionalization factor = (dU/dz) ^ 
time 

activation temperature 
temperature 

velocity difference across mixing layer 
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turbulence intensity 
velocity field 

free-stream speed of high-speed layer 

free-stream speed of low-speed layer 

streamwise, lateral, and spainwise coordinates 

corresponding velocity components 

initial concentration field scaling factor 

initial concentration field scaling factor 

mixing layer half-width 

mean product thickness 

circulation of the vortex 

mean vorticity thickness 

conserved scalar defined in Equation (12) 

longitudinal integral scale 

Taylor microscale 

kinematic viscosity 

density of fluid 

conserved scalar defined in Equation (16) 
stream function 

stream function for most unstable mixing layer mode 
stream function for subharmonic of most unstable mixing layer mode 
angular velocity 
characteristic length scale 


63 



ACKNOWLEDGEMENT S 


This work was supported by the NASA Lewis Research Center under Contract 
No. NAS3-23531. The code used in the 3D simulations was developed in 
collaboration with Dr. Steven A. Orszag. The work benefited from discussions 
with Professor R. E. Breidenthal and Dr. W.-R. Jou. 


64 



REFERENCES 


Acton, E. 1976. "The Modeling of Large Eddies in a Two-Dimensional Shear 
Layer", J. Fluid Mech., Vol . 76, Part 3, pp. 561-592. 

Ashurst, W. T., and P. K. Barr. 1981. "Lagrangian-Euler ian Calculation of 
Turbulent Diffusion Flame Propagation", presented at the Third Symp. Turb . 
Shear Flows, pp. 3.44-3.49, University of California, Davis. 

Batt, R. G. 1977. "Turbulent Mixing of Passive and Chemically Reacting Species 
in a Low-Speed Shear Layer", J. Fluid Mech. , Vol. 82, Part 1, pp. 53-95. 

Bilger, R. W. 1980. "Turbulent Flows with NonPremixed Reactants", in Turbulent 
Reacting Flows , edited by P. A. Libby and F. A. Williams, Springer-Ver lag. 

Borghi, R. 1974. "Chemical Reactions Calculations in Turbulent Flow", Adv . 

Geo phys ics , Vol. i8B, pp. 349-365. 

Breidenthal, R. E., Jr. 1978. "A Chemically Reacting, Turbulent Shear Layer", 
Ph.D. Thesis, California Institute of Technology. 

Broadwell, J. E., and R. E. Breidenthal, Jr. 1982. "A Simple Model of Mixing 
and Chemical Reaction in a Turbulent Shear Layer", J. Fluid Mech. , Vol 125, 
pp. 397-410. 

Browand, F. K. , and C.-M. Ho. 1983. "The Mixing Layer: An Example of Quasi 
Two-Dimensional Turbulence", to appear in J. de Mechanique . 

Browand, F. K. , and T. R. Troutt. 1980. "A Note on Spanwise Structure in the 
Two-Dimensional Mixing Layer", J. Fluid Mech. , Vol. 97, pp. 771-781. 

Brown, G. L., and A. Roshko. 1974. "On Density Effects and Large Structure in 
Turbulent Mixing Layers", J . Fluid Mech. , Vol. 64, pp. 775-816. 

Burke, S. P., and T. E. Schumann. 1928. "Diffusion Flames", Ind. Eng. Chem. , 
Vol. 120, p. 998. 

Chandrsuda, C., R. D. Mehta, A. D. Weir, and P. Bradshaw. 1978. "Effects of 
Free-Stream Turbulence on Large Structure in Turbulent Mixing Layers", J . 
Fluid Mech. , Vol. 85, pp. 693-704. 

Chorin, A. J. 1973. "Numerical Studies in Slightly Viscous Flow", J. Fluid 
Mech. , Vol. 57, pp. 785-796. 

Choudbury, P. R., M. Gerstein, and J. S. Cho. 1983. "A Novel Technique of 
Concentration Measurement in a Turbulent Reacting Medium", unpublished 
manuscript. 

Cullen, H. J. P. 1976. Quart. J. Roy. Meteorol. Soc. , 102, p. 77. 

Curl, R. L. 1963. "Dispersed Phase Mixing: 1. Theory and Effects in Simple 
Reactors", AIChE Journal , Vol. 9, pp. 175-181. 


65 



Danckwerts, P. V. 1952. Appl . Sci. Res. . Vol. 3, pp. 279-296. 


Deardorff, J. W. 1970. "A Numerical Study of Three-Dimensional Turbulent 
Channel Flow at Large Reynolds Numbers",. J. Fluid Mech., Vol. 41, 
pp. 453-480 

Deardorff, J. W. 1974. "Three-Dimensional Numerical Study of the Height and 
Mean Structure of a Heated Planetary Boundary Layer", Bound. Layer Met. , 

Vol. 7, pp. 81-106. 

Dimotakis, P., and G. L. Brown. 1976. "The Mixing Layer at High Reynolds 

Number: Large-Structure Dynamics and Entrainment", J. Fluid Mech. , Vol. 78, 
pp. 535-560. 

Donaldson, C. duP., and G. R. Hilst. 1972. "Effect of Inhomogeneous Mixing on 
Atmospheric Photochemical Reactions", Env. Sci. Tech., Vol. 6, pp. 812-816. 

Donaldson, C. duP., and A. K. Varma. 1976. "Remarks on the Construction of a 
Second-Order Closure Description of Turbulent Reacting Flows", Comb. Sci» 

Te ch . , Vol. 13, pp. 55-78. 

Dopazo, C. 1976. "A Probabilistic Approach to Turbulent Flame Theory", Acta 
Astronautica , Vol. 3, pp. 853-878. 

Forester, C. K. 1977. J. Comp. Phys . , Vol. 23, p. 1. 

Ghoniem, A. F., A. J. Chorin, and A. K. Oppenheim. 1982. "Numerical Modelling 
of Turbulent Flow in a Combustion Tunnel", Phil. Trans. R. Soc. Lond., Ser . 
A, Vol. 304, pp. 303-325. ” ‘ ‘ ~ “ 

Gibson, C. H. , and P. A. Libby. 1972. Comb. Sci. Tech. , Vol. 6, pp. 29-35. 

Haidvogel, D. B. , A. R. Robinson, and E. D. Schulraan. 1980. "The Accuracy, 
Efficiency, and Stability of Three Numerical Models with Application to 
Open Ocean Problems", J. Comp. Phys. , Vol. 34, pp. 1-34. 

Herring, J. R. 1974. "Approach of Axisymmetric Turbulence to Isotropy", Phys . 
FI. , Vol. 30, pp. 859-872. 

Herring, J. R., J. J. Riley, G. S. Patterson, Jr., and R. H. Kraichnan. 1973. 
"Growth of Uncertainty in Decaying Isotropic Turbulence", J ♦ Atmos . Sci . , 
Vol. 30, pp. 303-325. 

Hill, J. C. 1979. "Simulation of Chemical Reaction in a Turbulent Flow", in 
Proc. Second B. F. Ruth Chem. Eng. Res. Symp. , pp. 27-53, Iowa State 
University . 

Ho, C.-M., and L.-S. Huang. 1982. "Subharmonics and Vortex Merging in Mixing 
Layers", J . Fluid Mech . , Vol. 119, pp. 443-473. 

Hussain, A. K. M. F., and M. F. Zedan. 1978. "Effects of the Initial Condition 
on the Axisymmetric Free Shear Layer: Effects of the Initial Momentum 
Thickness", Phys . FI . , Vol. 21, pp. 1100-1112. 


66 



Janicka, J., and W. Kollman. 1979. "Prediction Model for the PDF of Turbulent 
Temperature Fluctuations in a Heated Round Jet”, presented at the 2nd Symp. 
Turb. Shear Flows, Imperial College, London, pp. 1.7-1.11. 

Jones, W. P., and J. H. Whitelaw. 1982. "Calculation Methods for Reacting 
Turbulent Flows: A Review", Combustion and Flame , Vol. 48, pp. 1~26. 

Karagozian, A. R. 1982. "An Analytical Study of Diffusion Flames in Vortex 
Structures", Ph.D. Thesis, California Institute of Technology. 

Konrad, J. H. 1977. "An Experimental Investigation of Mixing in Two- 

Dimensional Turbulent Shear Flows with Applications to Diffusion-Limited 
Chemical Reactions", Ph. D. . Thesis , California Institute of Technology. 

Libby, P. A., and F. A. Williams. 1981. "Some Implications of Recent 
Theoretical Studies in Turbulent Combustion", AIAA Journal , Vol. 19, 
pp. 262-274. 

Mans our , N. N., P. Moin, W. C. Reynolds, and J. H. Ferziger. 1979. "Improved 
Methods for Large-Eddy Simulations of Turbulence", in Turbulent Shear 
Flows I , p. 386, Springer-Verlag. 

Marble, F. E. 1982. "Growth of a Diffusion Flame in the Field of a Vortex", to 
appear in Luigi Crocco Anniversary Volume. 

Marble, F. E., and J. E. Broadwell. 1977. "The Coherent Flame Model for 
Turbulent Chemical Reactions", TRW Report 29314-6001-RU-00. 

Mason, H. B., and D. B. Spaulding. 1973. "Prediction of Reaction Rates in 
Turbulent Premixed Boundary Layer Flows", in Proc. First Symp. (European) 
on Comb . , pp. 601-606, Academic Press. 

McRae, G. J., W. R. Goodin, and J. : H. Seinfeld. 1982. "Numerical Solution of 
the Atmospheric Diffusion Equation for Chemically Reacting Flows", J. Comp. 
Phys . , Vol. 45, pp. 1-42. 

Mellor, A. M. , and C. R. Ferguson. 1980. "Practical Problems in Turbulent 
Reacting Flows", in Turbulent Reacting Flows , editied by P. A. Libby and 
F. A. Williams, Springer-Verlag. 

Metcalfe, R. W., and S. A. Orszag. 1974. "Numerical Simulation of Turbulent 
Jet Noise — Part 1", Flow Research Report No. 53. 

*Metcalfe, R. W., and J. J. Riley. 1981. "Direct Numerical Simulations of 
Turbulent Shear Flows", in Proc. Seventh Conf. Num. Methods in FI. Dyn. , 
pp. 279-284, Springer-Verlag. 

Michalke, A. 1964., "On the Inviscid Instability of the Hyperbolic-Tangent 
Velocity Profile", J. Fluid Mech. , Vol. 19, pp. 543-556. 

Moin, P., and J. Kim. 1982. "Numerical Investigation of Turbulent Channel 
Flow", J. Fluid Mech., Vol. 118, pp. 341-377. 


67 



Mularz, E. J. 1983. "New Trends in Combustion Research for Gas Turbine 

Engines", presented at the 6th Intern. Symp. Air Breath. Engines; also NASA 
T 83338; also AVRADCOM TR 83-C-l. 

Mungal , M. G. 1983. "Experiments on Mixing and Combustion with Low'Heat 

Release in a Turbulent Shear Flow", Ph.D. Thesis, California Institute of 
Technology. 

Mungal, M. G. , P. E. Dimotakis, and J. E. Broadwell. 1983. "Turbulent Mixing 
and Combustion in a Reacting Shear Layer", AIAA Paper No. 83-0473. 

Nielsen, K. A., and J. C. Hill. 1977.' "Numerical Simulation of Turbulent 
Mixing with Chemical Reaction", 70th Annual AIChE Mtg. 

Norton, 0. P. 1983. "The Effects of a Vortex Field on Flames with Finite 
Reaction Rates", Ph.D. Thesis, California Institute of Technology. 

O'Brien, E. E. 1969. "Postulate of Statistical Independence of Decaying 
Reactants in Homogeneous Turbulence", Phys. FI., Vol. 12, pp. 1999-2005. 

O'Brien, E. E. 1971. "On Turbulent Mixing of Two Rapidly Reacting Chemical 
Species", Phys . FI , Vol. 14, pp. 1326-1331. 

O'Brien, E. E. 1981. "Statistical Methods in Reacting Turbulent Flows", AIAA 
Journal , Vol. 19, pp. 366-371. 

Oran, E. S., and J. P. Boris. 1981. "Detailed Modeling of Combustion Systems", 
Prog. Energy Combust. Sci. , Vol. 7, pp. 1-71. 

Orszag, S. A. 1971. "Numerical Simulation of Incompressible Flows within 
Simple Boundaries: Accuracy", J. Fluid Mech. , Vol. 49, pp. 75-112. 

Orszag, S. A., and Y.-H. Pao. 1974. "Numerical Computation of Turbulent Shear 
Flows", Adv. Geophysics , Vol. 18A, pp. 225-236. 

Orszag, S. A., and G. S. Patterson, Jr. 1972. "Numerical Simulation of Turbu- 
lence", in Statistical Models and Turbulence , p. 127, Spr inger-Verlag. 

Patera, A. T., and S. A. Orszag. 1981. "Transition and Turbulence in Planar 
Channel Flows", in Proc. Seventh Conf. Num. Methods in FI. Dyn. , 
pp. 329-335, Springer-Ver lag. 

Peyret, P., and T. D. Taylor. 1983. Computational Methods for Fluid Flows , 

Spr inger-Verlag. 

Pope, S. B. 1979. Phil. Trans. Roy. Soc. Lond. , Ser. A, Vol. 291, pp. 529. 

Pratt, D. T. 1976. "Mixing and Chemical Reaction in Continuous Combustion", 
Prog. Energy Comb. Sci. , Vol. 1, pp. 73-86. 

Raymond, W. H. , and Gardner, A. 1976. Monthly Weather Rev. , 104, p. 1583. 


68 



Riley, J. J., and R. W. Metcalfe. 1978. "The Direct Numerical Simulations of 
the Turbulent Wakes of Axisymmetric Bodies — An Interim Report", Flow 
Research Report No. 135, NASA-CR-152282. 

*Riley, J. J., and R. W. Metcalfe. 1980a. "Direct Numerical Simulations of a 
Perturbed, Turbulent Mixing Layer", AIAA Paper No. 80-0274. 

*Riley, J. J., and R. W. Metcalfe. 1980b. "Direct Numerical Simulations of the 
Turbulent Wake of an Axisymmetric Body", in Turbulent Shear Flows II , 
pp. 78-93, Springer-Ver lag. 

Riley, J. J., and G. S. Patterson, Jr. 1974. "Diffusion Experiments with 
Numerically Integrated Isotropic Turbulence", Phys . FI , , Vol. 17, 
pp. 292-297 

Schumann, V., and G. S. Patterson, Jr. 1978. "Numerical Study of the Return of 
Axisymmetric Turbulence to Isotropy", J. Fluid Mech., Vol. 88, pp. 711-735. 

Spielman, L. A., and 0. Levenspiel. 1965. "A Monte-Carlo Treatment for 

Reacting and Coalescing Dispersed Phase Systems", Chem. Eng. Sci. , Vol. 20, 
pp. 247-254. 

Storch, H. 1978. Beitr'Age Phys. Atmosphare , 51, p. 189. 

Sturgess, G. J. 1983. "Aerothermal Modeling — Phase I", United Technologies 
Report CR 168202. 

Sutton, G. W. 1968. AIAA J. „ Vol. 6, pp. 1403-1405. 

Tennekes, H. , and J. L. Lumley. 1972. A First Course in Turbulence . MIT Press. 

Toor, H. L. 1962. "Mass Transfer in Dilute Turbulent and Nonturbulent Systems 
with Rapid Irreversible Reactions and Equal Dif fusivities" , AIChE Journal , 
Vol. 8, pp. 70-78. 

Toor, H. L. 1969. "Turbulent Mixing of Two Species With and Without Chemical 
Reactions", Ind. Eng. Chem. Fundam. , Vol. 8, pp. 655-659. 

Wallace, A. K. 1981. "Experimental Investigation of the Effects of Chemical 
Heat Release in the Reacting Turbulent Plane Shear Layer", Ph.D. Thesis, 

The University of Adelaide. 

Wax, N., Editor. 1954. Selected Papers on Noise and Stochastic Processes . 

Dover Publications, Inc. 

Wengle, H. , and Seinfeld, J. 1978. "Pseudospectral Solution of Atmospheric 
Diffusion Problems", J. Comp. Phys., Vol. 26, pp. 87-106. 

Williams, P. A. 1965. Combustion Theory . Addison-Wesley . 

Wygnanski, I., and H. E. Fiedler. 1970. "The Two-Dimensional Mixing Region", 

J. Fluid Mech. , Vol. 41, pp. 327-361. 


69 



Wygnanski, I., D. Oster, and H. Fiedler. 1980. "A Forced, Plane, Turbulent 
Mixing-Layer; A Challenge for the Predictor", in Turbulent Shear Flows 11 , 
pp. 314-326, Springer-Ver lag. 

Zalesak, S. T. 1979. "Fully Multidimensional Flux-Corrected Transport 
Algorithms for Fluids", J. Comp. Phys., Vol. 31, pp. 335-362. 

*Reprints available from the authors c/o Flow Research Company, Kent, 
Washington 98032. 


70 



Run No . 

e H 

TABLE 1. SIMULATION PARAMETERS* 

Two-dimensional simulations 
e SH ° E 

V 

Box 

MR01A 

0.1 

0 

1.00 

0.100 

1.00 

0.003 

2 tt 

MR01B 

o. 

0 

1.00 

0 . 100 

1.00 

0.003 

2 tt 

MR02A 

0.1 

0 

0.25 

0.100 

1.00 

0.003 

2 -n 

MR02B 

0 

0 

0.25 

0.100 

1.00 

0.003 

2 TT 

MR03A 

0.1 

0 

0.25 

0.200 

1.00 

0.003 

2 IT 

MR03B 

0 

0 

0.25 

0.020 

1.00 

0.003 

2 IT 

MR03C 

0.1 

0 

0.25 

0.010 

1.00 

0.003 

2 IT 

MR03D 

0 

0 

0.25 

0.010 

1.00 

0.003 

2 TT 

MR04A 

0.1 

0 

0.25 

0.005 

1.00 

0.003 

2 IT 

MR04B 

0 

0 

0.25 

0.005 

1.00 

0.003 

2 7T 

MR04C 

0.1 

0 

0.25 

0.020 

2.00 

0.003 

2 IT 

MR04D 

0 

0 

0.25 

0.020 

2.00 

0.003 

2 TT 

MY12B 

0.1 

0 

0.25 

0.005 

0 

0.003 

2 IT 

MY25A 

0.1 

0 

0.25 

0.005 

1.00 

0.003 

4 TT 

MY25B 

0.1 

0.1 

0.25 

0.005 

1.00 

0.003 

4 TT 

MY25C 

0 

0.1 

0.25 

0.005 

1.00 

0.003 

4 TT 

AG10A 

0.1 

0 

0.25 

0.005 

oo 

0.003 

2 IT 



Three- 

■dimensional simulations 



Run No . 

Realization 

z 

o 

D 

R 

V 

Box 

JN20B 

I 


0.25 

0.005 

1.00 

0.003 

2 TT 

JN24A 

I 


0.25 

0.01 

1.00 

0.003 

2 TT 

JL18A 

I 


0.25 

0.005 

1.00 

0.003 

4 TT 

JL26A 

II 

0.25 

0.005 

1.00 

0.003 

4 TT 

AG13A 

I 


0.25 

0.005 

0 

0.003 

2 IT 

AG13B 

I 


0.25 

0.005 

00 

0.003 

2 TT 

* The quantities in 

this table 

are nondimens ionalized by AU, the 

mean velocity 

difference across 

the mixing layer, 

T = (dU/dz) _i |t = 

0, the inverse mean 

velocity gradient at the centerline 

, and by the length scale L 

= UT. The 


columns ch and £SH specify whether the most unstable mode and/or its subhar- 
monic was introduced in the initial conditions. The value ejj = 0.1 produced 
an initial velocity field with w' = 0.071 U at the centerline. The initial 
concentration profile length scale is z 0 [see Equations (17) and (18)], D is 
the reactant diffusivity, R the reaction rate coefficient, v the kinematic 
viscosity, and Box is the length of a side of the computational domain, which 
was the same in all spatial directions. Re = Uz /v is the Reynolds number, 

Sc = n/D is the Schmidt number, D = z RC /U and°D = Rz 2 C / D are the 

.1 o . 00 II . o °o 

Darakohler numbers, DT is the time-step increment, and N is the number of 
Fourier modes in each spatial direction. Runs with ch = GSH = ® correspond to 
laminar, unidirectional (one-dimensional) flow. 
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TABLE 1. SIMULATION PARAMETERS (Cont.) 


Two-dimensional simulations 


Run No . 

Re 

Sc . 

D I 

D II 

DT 

N 

MR01A 

333 

0.030 

1.00 

10.000 

0.05 

64 

MR01B 

333 

0.030 

1.00 

10.000 

0.05 

64 

MR02A 

83 

0.030 

0.25 

0.625 

0.05 

64 

MR02B 

83 

0.030 

0.25 

0.625 

0.05 

64 

MR03A 

83 

0.015 

0.25 

0.313 

0.05 

64 

MR03B 

83 

0.150 

0.25 

3.125 

0.05 . 

64 

MR03C 

83 

0.300 

0.25 

6.250 

0.05 

64 

MR03D 

83 

0.300 

0.25 

6.250 

0.05 

64 

MR04A 

83 

0.600 

0.25 

12.500 

0.05 

64 

MR04B 

83 

0.600 

0.25 

12.500 

0.05 

64 

MR04C 

83 

0.150 

0.50 

6.250 

0.05 

64 

MR04D 

83 

0.150 

0.50 

6.250 

0.05 

64 

MY12B 

83 

0.600 

0 

0 

0.05 

64 

MY25A 

83 

0.600 

0.25 

12.500 

0.025 

128 

MY25B 

83 

0.600 

0.25 

12.500 

0.025 

128 

MY25C 

83 

0.600 

0.25 

12.500 

0.025 

128 

AG10A 

83 

0.600 

0 

0 

0.05 

64 




Three-dimensional simulations 



Run No • 

Re 

Sc 

°I 

D II 

DT 

N 

JN20B 

83 

0.600 

0.25 

12.500 

0.05 

64 

JN24A 

83 

0.600 

0.25 

12.500 

0.05 

64 

JL18A 

83 

0.600 

0.25 

12.500 

0.05 

64 

JL26A 

83 

0.600 

0.25 

12.500 

0.05 

64 

AG13A 

. 83 

0.600 

0 

0 

0.05 

64 

AG13B 

83 

0.600 

00 

CO 

0.05 

64 


1906R/1872R 
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Figure la t 


0 




Figure lc t = 16 

Figure 1. Plots of Lateral (y) Vorticity in an x-z Plane (Taken from Riley 
and Metcalfe, 1980a) 
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U(z) 
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Figure 4. Problem Geometry - Initial Conditions 
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Figure 6. Computed Concentration Field After One Revolution - No Diffusion 
or Reaction 
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Figure 7. Computed Concentration Field After One-Half Revolution - No 
Diffusion or Reaction 
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Figure 8. Computed Concentration Field After One-Half Revolution - No 
Diffusion, R « 1 .0 
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Figure 9. Computed Concentration Field After One-Half Revolution - No 
Diffusion, R 88 4.0 
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Figure 10. Constant Contour Plot of Initial Concentration Field Given by 
Equation (31) 
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Figure 11. Perspective Plot of Initial Concentration Field Given by 
Equation (31) 
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Figure 12. Constant Contour Plot of Computed Concentration Field After One 
Revolution - Diffusion but No Reaction 
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Figure 13. Perspective Plot of Computed Concentration Field After One 
Revolution - Diffusion but No Reaction 


Figure 14. Constant Contour Plot of Conserved Scalar Computed Indirectly 

from the Difference Between the Species Concentrations - t = 12 
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12 


Figure 15. Constant Contour Plot of Conserved Scalar Computed Directly 
from the Difference Between the Species Concentrations - t = 
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PRODUCT 



Figure 16. Total Product Versus Time for Two Different Values of the 
Reaction Rate Coefficient R - One-Dimensional Simulations. 
Also Shown Is the Infinite Reaction Rate Solution 
[Equation (37)] 


88 





DIFFUSIVITY 

Figure 17. Total Product Versus Diffusivity for a Fixed Time (t = 24) - 


One- and Two-Dimensional Single Rollup Simulations. Also Shown 
Is the Infinite Reaction Rate Solution [Equation (37)] 
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PRODUCT 





t 

Figure 18. Total Product Versus Time for Two Different Initial Concen- 
tration Profiles - One- and Two-Dimensional Single Rollup 
Simulations 
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Figure 19e t - 16 

Figure 19. Plots of Vorticity Contours for a Sequence of Times - Case 1 
(Fundamental Mode Alone) 
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Figure 20c t - 36 


lots of Vorticity Contours for a Sequence c 
Subharmonic Mode Alone) 
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Figure 22. Mean Velocity Half-Width, zj^ (- z^/2^’ Versus Time ~ 

Two-Dimensional Simulations (X Is the Wavelength of the Most 
Unstable Mode) 


Figure 23a Species A 



Figure 23b Species B 


Figure 23. Plots of Concentration Contours for Species A and B at 

t = 12 for Case 1 (Fundamental Mode Alone) - Two-Dimensional 
Simulations 
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Figure 24. Plot of Product Concentration Contours at t = 12 for Case 1 
(Fundamental Mode Alone) - Two-Dimensional Simulations 


96 



Figure 25. 


Plot of Product Concentration Contours at t = 24 for Case 1 
(Fundamental Mode Alone) - Two-Dimensional Simulations 
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Figure 27. 


Plot of Product Concentration Contours at t = 24 for Case 2 
(Subharmonic Mode Alone) - Two-Dimensional Simulations 
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Figure 28a t = 12 



Figure 28c t = 36 

Figure 28. Plots of Concentration Contours for Species A for a Sequence of 
Times for Case 3 (Fundamental and Subharmonic Added Together 
Out-of-Phase) - Two-Dimensional Simulations 
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Figure 29. Plot of Product Concentration Contours at t = 24 for Case 3 
(Fundamental and Subharmonic Added Together Out-of-Phase) 

- Two-Dimensional Simulations 
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Figure 30. Total Product Versus Time for Different Dif fusivities - One- 
and Two-Dimensional Simulations 
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REACTION RATE 



Figure 31. Total Reaction Rate Versus Time for Different Diffusivities 
One- and Two-Dimensional Simulations 
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MIXEDNESS PARAMETER 
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Figure 32. Mixedness Parameter Versus Time for Different Dif fusivities - 
Case 1 (Fundamental Mode Alone), Two-Dimensional Simulations 
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PRODUCT 



t 

Figure 33. Total Product Versus Time - Two-Dimensional Simulations 
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MIXEDNESS PARAMETER 
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Figure 34. Mixedness Parameter Versus Time - Two-Dimensional Simulations 
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PRODUCT THICKNESS 
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Figure 35. Nondimensional Product Thickness Versus Time - Two-Dimensional 
Simulations 


107 






IU 



Figure 37. RMS of the Concentration Fluctuation of Species A Versus z for 
Several Different Times for Case 1 (Fundamental Mode Alone) ~ 
Two-Dimensional Simulations 
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Figure 38. RMS of the Concentration Fluctuation and Average Concentration 
of Species A at t = 12 for Case 1 (Fundamental Mode Alone) - 
Two-Dimensional Simulations 
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Figure 39. Average Product Concentration Versus z for Several Different 
Times for Case 1 (Fundamental Mode Alone) - Two-Dimensional 
Simulations 
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Figure 40. The Average of the Product of the Concentrations, the Product 

of the Averages, and the Concentration Correlation Versus z for 
Case 1 (Fundamental Mode Alone) at t = 12 - Two-Dimensional 
Simulations 




Figure 41. Various Concentration Moments Versus z for Case 1 (Fundamental 
Mode Alone) at t = 12 - Two-Dimensional Simulations 






PROBABILITY 



Figure 42. The Probability Density Function for the Concentration of 

Species A at the Transverse Distance — it/ 2 for a Sequence of 
Times for Case 1 (Fundamental Mode Alone) - Two-Dimensional 
Simulations 


114 


PROBABILITY 
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The Probability Density Function for the Concentration of 
Species A at t = 12 for Various Transverse Distances for Case 1 
(Fundamental Mode Alone) - Two-Dimensional Simulations 




Figure 44a t = 6 


Figure 44b t = 12 



Figure 44c t = 18 Figure 44d t = 24 


Figure 44. Plots of Lateral (y) Vorticity Contours for a Sequence of Times 
in an x-z Plane - Single Rollup Case, Three-Dimensional 
Simulations, Run JN20B 


116 









Figure 45c t - 18 Figure 45d t - 24 


Figure 45. Plots of Lateral (y) Vorticity Contours for a Sequence of Times 
in an x-z Plane Located a Lateral Distance Ly/2 from that 
Depicted in Figure 44 - Single Rollup Case, Three-Dimensional 
Simulations, Run JN20B 
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MEAN VELOCITY HALF-WIDTH 
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Figure 47. Mean Velocity Half-Width, z^ (= Versus Time - 

Three-Dimensional Simulations 
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Figure 49. Mean Prc 
Simulat i 
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Figure 50. Nondimensional Product Thickness Versus Equivalence Ratio 
(Taken from Mungal, 1983) 


NONDIMENSIONAL PRODUCT THICKNESS 




t 

Figure 51. Nondiraensional Product Thickness Versus Time - 
Three-Dimensional Simulations 
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Figure 52. Mixedness Parameter Versus Time - Three-Dimensional Simulations 
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Figure 54. Average Concentration of Species B Versus z for a Sequence of 
Times - Realization I, Double Rollup Case, Three-Dimensional 
Simulations 
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Figure 56. RMS of the Fluctuating Concentration of Species A Versus z for 
a Sequence of Times - Realization I, Double Rollup Case, 
Three-Dimensional Simulations 
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Figure 58.' RMS of the Fluctuating Concentration of Species A Versus z for 
a Sequence of Times - Realization II, Double Rollup Case, 
Three-Dimensional Simulations 
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Figure 61 
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. Various Concentration Moments Versus z at t - 12 - Realization II, 
Double Rollup Case, Three-Dimensional Simulations 
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Figure 62. Similarity Plots of C A / Versus z/z M (= z/z^/2^ for Various 
Times - Double Rollup Case, Finite Reaction Rate, Three- 
Dimensional Simulations 
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Figure 63. Similarity Plots of Cg/C^ Versus z/z M ( = z / z l/2^ f° r Var i° us 
Times - Double Rollup Case, Finite Reaction Rate, Three- 
Dimensional Simulations 


135 





Figure 64. Similarity Plots of C A /C Aoo Versus z/z M (= z/z^) for Various 
Times - Double Rollup Case, Finite Reaction Rate, Three- 
Dimensional Simulations 
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Figure 71. RC^Cg/ [min(C^, Cg) /T] Versus z at t = 12 - Double Rollup Case, 
Finite Reaction Rate, Three-Dimensional Simulations 
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Figure 72. Min(C^, Cg) /C^Cg Versus Time Measured at z = 0 - Double Rollup 
Case, Finite Reaction Rate, Three-Dimensional Simulations 
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Figure 73. c a c B^A<» Ver sus z at t = 16 - Case 1 (Fundamental Mode Alone), 

Both Zero and Finite Reaction Rates, Two-Dimensional Simulations 
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